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ABSTRACT 

Two  algorithms  for  generating  a ncn-homogeneous 
Poisson  process  with  log-quadratic  intensity  function 
A (t)  = exp  (aQ  ♦ c^t  + a2t2)  S>are  implemented  into 
computer  programs  and  compared  for  relative  speed, 
core  storage  requirements  and  fidelity.  By  simulating 
several  cases  of  non-homogeneous  Poisscn  processes 
with  log-quadratic  intensity  functions  it  is  shewn 
that  the  Poisson- decomposition  and  gap  statistic 
algorithm, — developed  by  Professor  P.A.W.  Lewis,  Haval 
Postgraduate/ School,  Monterey,  California,  and  G.S. 
Shedler,  / IBM  Research  Laboratory,  San  Jose, 
Ca  lif  crniaVNsubstantiall  y reduces  computation  time 
from  that  required  by  an  algorithm  that  uses  a 
time-scale  transformation  of  a homogeneous  Poisson 
process.  -—-The  faster  algorithm  employs  a rejection 
technique  in  exjunction  with  a method  for  simulating 
the  nen-homogeneous  Poisson  process  with  intensity 
function  A(t)  = $xp  ( t)  by  generation  of  gap 
statistics.  Although  additional  core  storage  is 
required  by  thgf  Lewis  and  Shedler  algorithm,  the 
resulting  gaib  in  computing  efficiency  is  so 
significant  that  it  outweighs  the  memory 
consideration.  The  experience  gained  from 
implementing  the  algorithm  has  led  to  several 
possibilities  which  are  suggested  for  improving  the 
efficiency  of  the  Poisson-decoraposition  and  cap 
statistic  algorithm. 
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INTRODUCTION 


Many  faniliar  physical  and  operational  processes  are 
well  described,  in  whole  or  part,  by  examining  their  "event 
streams"  over  time.  Some  such  well-known  processes  are:  the 
flow  of  traffic  through  an  intersection;  the  arrival  of 
persons  at  seme  service  facility  such  as  a bank  teller's 
window,  a service  station  fuel  pump  or  a grocery  store  check 
out  counter;  and,  the  arrival  of  telephone  calls  or  radio 
transmissions  at  some  switchboard  or  other  type  of 
communications  terminal.  Analogous  processes  abound  both  in 
nature  and  in  the  course  of  our  everyday  lives.  Ey  the 
proper  definition  of  an  "event"  in  these  situations,  the 
process  being  observed  will  be  characterized  by  the 
probabilistic  nature  of  the  flow,  over  time,  of  the  events 
of  which  it  is  composed.  In  the  above  three  examples  the 
events  could  be  defined  respectively  as  the  arrival  of  a 
vehicle  at  the  intersection,  the  arrival  of  a customer  at 
the  service  facility,  and  the  arrival  of  the  telephone  call 
or  radio  transmission  at  the  terminal.  The  process  may  then 
be  analyzed  by  examining  the  interaction  of  various  event 
streams  with  different  intersection  configurations,  service 
policies  and  terminal  capacities. 

A common  method  used  to  perform  such  analyses  is  to 
consider  the  event  streams  to  be  homogeneous.  This  could 
mean  that  the  expected  number  of  events  to  occur  in  any  two 
or  more  time  intervals  of  equal  length  is  the  same  (simple 
homogeneity),  or  that  the  distribution  of  the  number  of 
events  occurring  in  any  two  or  more  equal  time  intervals  is 
the  same  (complete  homogeneity) . The  homogeneous  Poisson 
process  is  often  used  as  a tool  in  the  analysis  of  such 
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activities.  For  very  simple  systems  involving  event  streams 
the  use  of  the  homogeneous  Poisson  process  as  a model  leads 
to  tractable  analytical  results.  Most  systems  of  interest, 
however,  are  not  amenable  to  purely  analytical  methods,  and 
simulation  of  the  processes  by  digital  means  becomes 
necessary. 

The  assumption  of  homogeneity  in  event  streams  is  often 
a very  restrictive  one.  The  "rush  hour"  phenomenon, 
well-known  and  abundantly  cursed  by  motorists,  provides 
cogent  evidence  that  the  modeling  of  event  streams  is  not 
always  well  served  by  the  homogeneity  assumption.  The 
intensity  of  event  streams  varies  over  tine  for  many 
physical  cr  operational  processes.  In  these  cases,  purely 
analytical  methods  must  be  abandoned  almost  immediately  in 
favcr  of  simulation  techniques.  (The  intensity  of  the  event 
stream  is  defined  to  be  the  derivative  with  respect  to  t of 
the  expected  number  of  events  in  an  interval  of  length 

<o,  1 3. ) 

The  ncn-hcmogeneous  Poisscn  process  is  often  employed  in 
the  analysis  of  processes  that  exhibit  gross  departures  from 
the  homogeneous  event  stream  criterion.  If  these  processes 
are  to  be  simulated,  it  is  necessary  to  first  describe  the 
nature  of  the  inhomogeneity  (i.e.  how  does  the  intensity  of 
the  event  stream  vary  over  time?)  and  then  to  artificially 
generate  event  streams  that  behave  in  accordance  with  the 
descripticn. 


Of  course  there  are  infinite  variations 
inhomogeneities  that  can  occur  or  that  can 
However  it  is  intuitively  appealing  to 
streams  that  display  varying  intensities  of 
types : 


in  the  types  of 
be  construed, 
consider  event 
the  following 


i)  increasing  continuously  over  time; 
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ii)  decreasing  continuously  over  time; 

iii)  increasing  and  then  decreasing  continuously  over 
tine ; 

iv)  decreasing  then  increasing  continuously  over  time. 

& ncn-hcaogeneous  Poisson  process  that  has  quadratic 
properties  can  be  aanipulated  to  produce  the  above- aentioned 
effects.  The  effective  siaulation  of  such  a process  on  the 
computer  is  the  subject  of  this  thesis  and  has  been 
motivated  by  the  worfc  of  P.  A.  W.  Lewis,  Professor,  Naval 
Postgraduate  School,  and  G.  S.  Shedler,  IBM  Research 
Laboratory. 

There  are  of  course  other  types  of  inhomogeneities,  such 
as  cyclic  variations  (tiae  of  day  effect),  but  we  do  not 
consider  them  here.  They  have  been  discussed  by  COX 
[Ref.  2]  and  LEWIS  [Ref.  9].  In  particular  LEWIS  [Ref.  9] 
describes  a process  consisting  of  arrivals  at  an  intensive 
care  unit  in  a hospital.  It  is  shown  empirically  that,  in 
addition  to  the  time  of  day  cycle,  long  tern  fluctuations  in 
the  intensity  function  can  be  adequately  described  by  an 
intensity  function  whose  logarithm  is  quadratic. 

LEWIS  and  SHEDLER  [Ref.  11]  proposed  a new  method  for 
generating  a non-homogeneous  Poisson  process  with  an  event 
stream  intensity  (rate)  function  that  is  of  degree-two 
exponential  polynomial  form.  (The  use  of  exponential 
polynomials  is  natural  in  this  context  since  an  intensity 
function  is  a positive  function.)  The  new  method  appears  to 
have  the  virtue  of  increased  efficiency  over  the  more 
conventional  time-scale  transformation  technique  when 
implemented  on  a high  speed  digital  computer. 

Efficiency  in  this  context  is  measured  in  terms  of 
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conputer  memory  requirements  and  computational  speed.  The 
problem  of  efficiency  comparison  is  recognized  by  LEWIS  and 
SHEDLER  in  the  final  pages  of  their  paper  [Ref.  11,  p.  15]: 
"There  remains  the  question  of  efficiency  (of  the  proposed 
algorithm)  . . . for  generation  of  a non-homogeneous  Poisson 
process  with  degree-two  exponential  polynomial  rate 
function,  relative  to  generation  via  time  scale 
transformation  of  a homogeneous  Poisson  process." 

After  seme  brief  discussion  of  the  requirements  of 
implementing  the  time-scale  transformation  algorithm  the 
report  concludes  . . . "We  therefore  would  expect  the  exact 
method  of  (the  proposed  algorithm)  to  be  much  faster, 
although  at  the  expense  of  some  complexity  of  programming." 

The  objective  of  this  author  has  been:  to  implement 
both  the  algorithm  of  LEWIS  and  SHEDLER  [Ref.  11]  and  the 
conventional  time-scale  transformation  algorithm  on  the  IBM 
360/67  Computer  System  in  FORTRAN  IV  language;  to  define 
reasonable  measures  of  effectiveness  for  comparing  the  two 
algorithms;  and  to  determine  which  algorithm  is  the  more 
efficient  in  terns  of  the  measures  of  effectiveness  defined. 


Section  II  discusses  the  definition  of  a non-homogeneous 
Poisson  process.  It  also  states  some  special  properties  of 
Pcisscn  processes  that  are  used  in  the  development  of  the 
algorithms  investigated  in  this  thesis,  and  concludes  with  a 
general  discussion  of  two  basic  methods  of  generating  a 
non-homogenecus  Poisson  process.  Section  III  gives  a step 
by  step  description  of  the  two  algorithms  that  were 
implemented  into  computer  programs.  The  method  of 
implementation  is  the  topic  of  Section  17.  Methods  of 
comparing  the  algorithms  are  presented  in  Section  /. 
Section  VI  presents  conclusions  drawn  from  the  comparison 
and  makes  a recommendation  for  improving  the  LEWIS  and 
SHEDLER  [Ref.  11]  algorithm.  Other  recommendations  for 
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further  study  are  also  listed  in  Section  VI.  Appendixes 
provide  additional  details  that  would  have  been  awkward  to 
include  in  the  main  body  of  the  thesis.  Computer  program 
listings  are  provided  after  Appendix  E. 
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II.  DEFINITION  AND  PROPERTIES  OF  NON-HOKOGENEOUS  PCISSON 

PROCESSES 


A.  GENERAL 

This  section  will  present  basic  definitions  and 
explanations  concerning  the  concept  of  non-homogeneous 
Poisson  processes.  References  cited  nay  be  consulted  if  a 
more  in-depth  understanding  is  desired.  Only  the 
fundamental  concepts  necessary  for  understanding  the 
specific  non-homogeneous  Poisson  process  under  consideration 
are  presented. 

B.  DEFINITION  OF  A NON- HOMOGENEOUS  POISSON  PROCESS 

LEWIS  and  SHEDLER  [Ref.  11]  define  the  non-hooogeneous 
Poisson  process  on  the  real  line  as  follows: 

1.  The  number  of  events  in  any  finite  set  of 
non-overlapping  intervals  are  independent  random  variables. 

2.  Let  A (t)  be  a monotone  non- deer  easing 
right-continuous  function  which  is  bounded  on  any  finite 
interval.  Then  the  number  of  events  in  any  interval,  e.g. 
(0,tg),  has  a Poisson  distribution  with  parameter 

A(tQ)  - A (0)  . 

The  function  A(t)  - A (0)  is  called,  among  other  things,  the 
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mean  value  function  since  the  expected  number  of  events  in 
an  interval  {0,t]  is  just  A (t)  - A (0)  . Property  1 is  the 
independent  increment  property  of  the  Poisson  process;  it  is 
basic  to  the  idea  of  a Poisson  process.  Figure  1 provides  a 
graphical  representation  of  the  definition  above. 

The  above  definition  insures  that 

( A(t)  - A (0) } / { A (tQ)  - A (0) } meets  the  following 
criteria  for  an  arbitrary  function  F (t)  , to  be  a valid 
distribution  function  on  the  interval  (0,tg),  c.f.  LARSON 


[Ref. 

7,  ch.  3]; 

i) 

0 £ F ( t ) <_  1 

ii) 

lim  F (t ) = 0;  lim 

t 0 t ■*  tg 

F(t)  = 1, 

0 < 

t < 

iii) 

F (a)  <_  F (b)  for  all 

a <_  b in 

(0 

iv) 

lim  F (b  + h)  = F (b) 

for  all  b 

in 

(0 

h -*■  0 

where  h > 0 . 


Letting  F (t)  = { A (t)  - A (0)  } / [ A (tQ)  - A (0) } it  follows 
that  if  A (t)  is  absolutely  continuous  in  (0,tg]  then 
dF(t)/dt  = X(t)/{  A (tQ)  - A < 0) ) is  a valid  density  function 
on  the  interval  (0,t  ].  The  function  X (t)  is  called  the 
intensity  function  (or  rate  function)  of  the  process  and 


X(t)  = E{number  of  events  in  (0 , tQ  ] } 

= rjL  { A (t)  - A (0  ) } 
d 

= dt  A(t)  • 
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Situation:  Events  occur  randomly  In  time  (i.e.  along  t-axis) 

1)  If  X * number  of  events  In  fixed  Interval  (tj.t^J 

Y * number  of  events  In  fixed  Interval  (t^.tj] 

l * number  of  events  In  fixed  Interval  (tj.t^] 

S ■ number  of  events  In  fixed  Interval  (O.tpj; 

Then  X,  Y and  l must  be  Independent  random  variables  (note: 

S and  l are  not  Independent  because  of  overlap  In  the  Intervals). 

2)  Given  the  monotonlc  Increasing  right  continuous  function 
A ( t ) . the  number  of  events  In  any  fixed  Interval  (t^.tj] 

must  have  the  Poisson  distribution  with  parameter  A(t^)  - Aftj). 

X - Po1sson(A(tz)  - A(tj)) 

Y ~ PolSsoniA(tj)  - A(t2)} 
l - P*1»son(A(t4)  - a( tj) ) 

S ~ Pol sson{A( tp)  - A (0 ) > . 

Figure  1 - DEFINITION  OF  A NON-HONOGEN  SOUS  POISSON  PROCESS 

(GRAPHICAL  REPRESENTATION) 
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C.  THE  INTENSITY  PONCTION 


The  nest  intuitively  appealing  way  to  think  about  a 
non-homogeneous  Poisson  process  is  to  consider  the  variation 
in  the  intensity  of  the  event  stream  over  time.  The  concept 
of  an  intensity  function  whose  integral  meets  the  criteria 
of  the  definition  in  paragragh  B above  is  essential  tc  the 
modeling  and  simulation  of  a non-homogeneous  Poisson 
process.  When  one  starts  with  the  existence  of  1 (t) , the 
A(t)  is  often  called  the  integrated  intensity  (or  rate) 
function. 

The  intensity  functuion  reveals  the  instantaneous  rate 
of  arrivals  (in  the  event  stream)  as  a function  of  time. 
(This  function  must  be  a positive  function.)  For  example, 
if  telephone  calls  arrive  at  a switchboard  at  the  rate  of  5 
per  hour  at  0900  (9:00  a.m.),  increase  to  a peak  rate  of  20 
per  hour  at  1300  (1:00  p.m.),  then  decrease  to  a rate  of  5 
per  hour  at  1700  (5:00  p.m.),  the  intensity  function  could 
lock  something  like  Figure  2a.  Then,  by  plotting  the 
integral  of  the  intensity  function  over  the  interval  of 
interest  (i.e.  from  0900  to  1700)  it  is  obvious  that  a 
mcnotone-increasing,  right-continuous  and  bounded  function 
is  obtained  (see  Figure  2b)  . If  the  assumption  is  made  that 
the  arrival  stream  is  a Poisson  process,  i.e.  has 
independent  increments,  then  the  number  of  calls  received  in 
any  chosen  interval  (e.g.  0900-1000,  1230-1315,  1107-1632) 
is  distributed  as  a Poisson  random  variable  with  parameter 
equal  tc  the  difference  between  the  integral  evaluated  at 
the  right  end  point  of  the  interval  and  the  integral 
evaluated  at  the  left  end  point  of  the  interval.  These 
values  may  be  read  directly  from  Figure  2b.  Specifically, 
the  number  of  calls  received  in  an  eight-hour  working  day  is 
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a Fcissor  random  variable  with  parameter  120. 

Although  it  is  unlikely  that  telephone  calls  would 
arrive  at  a switchboard  in  accordance  with  the  convenient 
parabolic  intensity  function  of  Pigure  2a,  the  example 
serves  to  illustrate  two  important  points.  Pirst,  it  would 
not  be  realistic  to  assume  that  the  arrival  rate  of 
telephone  calls  at  a switchboard  would  be  constant 
throughout  a working  day.  Some  function  that  describes  an 
initially  increasing  and  finally  decreasing  rate  of  arrivals 
seems  more  akin  to  reality.  The  importance  of  being  able  to 
model  a non-homo geneous  Poisson  process  is  therefore 
established. 

Secondly,  it  is  obvious  that  the  definition  of  a 
ncn-homogeneous  Poisson  process  is  not  always  used  as  a 
starting  point  for  modeling  operational  processes.  Event 
streams  are  usually  thought  of  in  terms  of  their  underlying 
intensity  functions.  The  idea  of  an  intensity  function 
applies  to  any  model  for  an  event  stream,  and  is  not 
specific  to  a Pcisson  process.  The  further  step  of  modeling 
the  physical  process  as  a non-homogeneous  Poisson  process  by 
assuming  that  the  process  has  independent  increments  is 
taken  either  on  the  basis  of  empirical  evidence  or  physical 
reasoning.  Testing  for  independent  increments  in  a point 
process  is  discussed  in  COX  and  LEWIS  [Ref.  3,  ch. .6]  and 
LEWIS  [Ref.  9].  The  main  physical  reason  for  assuming 
independent  increments,  and  therefore  a Poisson  process,  is 
that  the  operational  process  is  the  superposition  of  many 
individual  event  streams.  For  instance,  in  a computer 
center  eguipped  with  several  interactive  time-sharing 
terminals,  the  event  stream  of  users  at  each  specific 
terxinal  night  be  assumed  to  be  an  arbitrary  point  process 
with  a certain  intensity  function.  The  event  stream  seen  by 
the  central  processing  unit  is  then  the  sum  total  (or 
superposition)  of  the  event  streams  of  the  individual 
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d,  if  there  are  enough  terminals,  should  be 
y Poisson.  This  property  of  superposition  of 
notions  is  discussed  in  the  next  paragraph. 


time  of  day 
11s  at  a Switchboard 


NOTIONS  FOR  TELEPHONE 


D.  DECOMPOSITION  AND  SUPERPOSITION  OF  POISSON  PROCESSES. 


To  obtain  a Poisson  process  with  intensity  function 
X (t)  = A^ft)  + X2  (t)  , we  may  superpose  two  Pcisson 
processes,  each  of  intensity  X^(t)  and  X2('t). 

We  may  decompose  the  intensity  function  of  any  event 
stream  into  two  or  more  component  event  streams.  However  if 
we  then  superpose  the  component  streams,  we  will  recover  the 
original  type  process  only  if  we  started  with  a Pcisson 
process.  For  example  begin  with  a renewal  process  with 
intensity  v (t)  = and  let  v =v1  ♦ v2  . If  two  renewal 
processes  with  intensities  and  are  superposed,  the 
resulting  process  is  not  a renewal  process. 

This  ucigue  property  may  be  exploited  when  simulating 
Poisson  processes.  It  permits  the  partitioning  cf  the 
intensity  function  to  take  advantage  of  any  special 
properties  cf  its  component  parts.  This  is  the  basis  for 
the  method  used  in  the  poisson-Decomposition  and  Gap 
Statistic  Algorithm  discussed  in  later  sections. 


E.  TWO  EASIC  METHODS  CF  GENERATING  A NON- HOMOGENSO US 
PCISSON  FBOCESS 


1 . Time-Scale  Transformation  Algorithm 


Consider  the  non-homogeneous  Poisson  process  with 
intensity  function  (t)  > 0,  0<t<tQ,  on  the  interval 

(0,tg].  The  integral  of  the  intensity  function  is  then 
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A (t)  and  the  number  of  events  in  the  interval  is  a Poisson 

random  variable  with  parameter  A (t  ) - A (0)  = (or 

A ( tQ ) * since  A (0)  = /®X(t)dt  = 0).  Now  if  T* , T* 

are  events  in  a unit  homogeneous  Poisson  process  (i.e.  a 

Poisson  process  with  constant  intensity  function 

X • (t)  = X'  = 1)  then  A-1  (T*) , ...,A-1(T*)  are  events  in  the 

1 n 

non-homogenecus  Poisson  process. 

This  result  gives  a procedure  for  simulating  a 

ncn-homogeneous  Poisson  process,  starting  with  a homogeneous 

Poisson  process,  which  is  analogous  to  the  probability 

integral  transform  method  of  producing  random  samples  from  a 

continuous  distribution  with  distribution  function  F(x)  when 

starting  with  uniformly  distributed  random  variables.  The 

latter  method  is  essentially  that  if  the  inverse  of  F(x)  can 

be  found,  then  by  generating  uniform  (0,1)  variates  u^ , 

. . . , u the  values 
n 

X.  — F (u.),  . . . / x — F (u  ) 

11  n n 

comprise  a sample  of  variates  from  the  desired  distribution. 

The  right-continuous  monotone  increasing  function 

A (t)  describing  the  ncn-homogeneous  Poisson  process  on  the 

interval  (0,to]  can  be  thought  of  as  a distribution  function 

which  has  teen  "scaled"  by  the  factor  u q • (Since 

A<V  = yQ/  then  A ) /y  * 1#  thus  A(t)/UQ  is  a valid 

distribution  function  on  (0,t  ].) 

0 

Tc  implement  the  time-scale  transformation  procedure 
one  can  use  the  following  basic  result  for  homogeneous 
Pcisson  processes:  given  that  n events  have  occurred  in  a 
homogeneous  Pcisson  process  over  a fixed  interval  ( C , t q] , 
those  events  are  uniformly  distributed  on  the  interval.  A 
prcof  of  this  property  is  given  in  PARZEN  [Ref.  14,  p.  140]. 
Therefore,  if  events  are  generated  as  a unit  Poisson  process 
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on  an  interval  of  length  Pg  by  first  obtaining  a 
Poisson  (Pg)  random  variate  n,  and  then  letting  the  order 
statistics  from  a random  sample  of  uniform  (0,Pg)  variates 
be  the  times  to  events  in  the  unit  Poisson  process,  and  the 
times  to  these  events  are  then  transformed  b^  the  inverse  of 
the  integrated  intensity  function,  i.e.  A (•) , the  effect 
is  the  same  as  that  obtained  from  the  probability  integral 
transform  method.  Pigure  3 illustrates  this  method 
graphically  and  also  provides  a flow  chart.  Note  the 
difference  however  between  this  procedure  and  the 
probability  integral  transform  procedure.  In  generating  a 
unit  Poisscn  process  by  this  method  we  need  a sample  whose 
size  is  random  (and  could  be  zero)  i.e.  Poisson  ( Pq)  . The 
probability  integral  transform  method  simply  transforms  a 
fixed  number  of  uniform  (0,1)  variates  into  variates  from 
seme  other  distribution. 


In  Appendix  A it  is  proved  that  scaling  both  the 
distribution  function  F(x)  and  the  interval  (0,1)  by  the 
same  factor  dees  not  affect  the  validity  of  the  probability 
integral  transform  method. 


The  unit  homogeneous  Poisson  process  may  also  be 
generated  by  adding  a sequence  of  unit  exponential  variates 
(variates  from  an  exponentially  distributed  random  variable 
with  mean  = 1)  until  the  sequence  of  partial  sums  of  the 
random  variables  first  exceeds  Pg  (see  Appendix  D) . This 
accomplishes  two  things.  First  it  provides  an  ordered 
sequence  of  events  from  a homogeneous  Poisson  process  on  the 
interval  (0,Pg].  Second,  it  determines  a realization  of  the 
Poisson  random  variable  N t with  parameter  A(tg)  = Pg  . Note 
that  given  n,  the  times  to  events  are  uniform  order 


the 


Figure  3 - GRAPHICAL  REPRESENT ATION  OF  TIME- 
TRANSFORMATION  METHOD 


The  difference  between  these  methods  for  generating 
a homogeneous  Poisson  process  should  also  be  noted.  The 
first  method  requires  a Poisson  variate  and  ordering  of  that 
number  of  uniform  random  variables.  The  second  requires 
generation  of  independent  exponential  variates.  The  second 
method  is  probably  most  basic  in  that  it  requires  only 
exponentially  distributed  random  variables,  and  these  can  be 
obtained  from  uniform  variates  by  the  inverse  probability 
integral  transform,  which  is  just  a logarithm.  The  method 

is,  however,  not  always  the  most  efficient. 

/ 

2 • Conditioning  and  ' Orde r Statistics  Algorithm  for  a 
Poisson  Process 


This  method  requires  the  result  of  a theorem 
sketched  by  LEWIS  and  COX  [Ref.  3,  ch.  2]  and  restated  here 
for  convenience  and  continuity;  it  is  an  extension  of  the 
result  on  conditioning  in  a homogeneous  Poisson  process 
which  results  in  a conditional  uniform  distribution  of  the 
times  to  events. 


Theorem  _1 : Assume  that  a non-homogeneous  Poisson  process  is 
observed  for  a fixed  time  (0,tg],  so  that  the  number  of 
events  in  (0,tg],  ^tg'  ^as  a P°isson  distribution  with 


parameter  A(tQ)  - A (0)  = Uq 


Then  if  T 


1' 


,T  denote 
n 


times-to-events  for  the  process  in  (0,tQ],  and  if  = n, 
conditional  cn  having  observed  n (>0)  events  in  (0,tQ],  the 
T^’s  are  distributed  as  order  statistics  from  a sample  with 
distribution  function 


F(t) 


. A ( t ) - A (0) 
A ( 1 0 ) — A ( 0V 


This  theorem 


reduces  the  problem  of  simulating  a 


non-homogeneous  Poisson  process  to  that  of  generating  a 
Pcisson  number  of  order  statistics  from  a fixed  distribution 
function.  That  is,  given  an  intensity  function  over  an 
interval  (a,b],  whose  definite  integral  A (t)  * /tX(s)ds  is 

di 

bounded  and  right-continuous  on  the  interval,  an  ordered 
sample,  fxcm  a population  with  distribution  function 


A (t)  - A (a) 
A (b ) - A (a) 


a < t < b (1) 


will  yield  the  desired  non-homogeneous  Poisson  process 
defined  by  the  intensity  function  X(t)  on  (a,b].  For 
simplicity  the  interval  will  hereafter  be  assumed  to  have 
its  left  end  point  at  zero  (a  = 0)  and  its  right  end  point 
at  some  arbitrary,  but  fixed  point  tQ,  (b  * tQ) . Using  this 
(0 , tg  ] interval  results  in  no  loss  of  generality,  and  (1) 
becomes  identical  with  the  expression  in  the  theorem. 


Many  methods  exist  for  obtaining  the  necessary  order 
statistics.  The  inverse  integral  transform  explained  above, 
decomposition  of  the  density  function  (see  LEWIS  Bef . 8), 
or  the  rejection-acceptance  method  discussed  later  in 
Section  IV-D,  are  all  possibilities. 

For  the  family  of  intensity  functions  addressed  in 
this  thesis  the  non-hcmogeneous  Poisson  process  is  obtained 
by  a combination  of  Poisson  decomposition,  an  algorithm  of 
LEWIS  and  SEEDLE2  [Ref.  13]  for  obtaining  a non-homogeneous 
Pcisson  process  with  log-linear  intensity  function,  and  the 
conditioning  and  order  statistics  theorem  given  above.  The 
procedure  involves  four  basic  steps. 

1.  The  intensity  function  is  decomposed  into  two 
ccmpcnents;  i.e.  X(t)  = _X(t)  ♦ X*  (t)  . 
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2.  A gap  statistics  algorithn  is 
a ncn-ho aogeneous  Poisson  process 
components,  A(t),  which  is  chosen 
a special  structure  (log-linear) . 


used  to  generate 
from  one  of  the 
such  that  it  has 


3.  A rejection-acceptance  routine  is  used  to 
produce  a sample  from  the  remaining  component, 
X*  (t) , i.e.  a Poisson  number  of  ordered  variates 
are  generated.  This  algorithm  is  described  in 
detail  in  Section  III-C.  This  sample,  when 
ordered,  becomes  a Poisson  process  also. 

4.  The  Poisson  processes  of  the  component 
intensity  functions  are  then  superposed  to  produce 
the  desired  non-ho mogeneous  Poisson  process. 

Figures  4,  5 and-  6 (Section  III)  illustrate  the  four  general 
steps  of  this  procedure. 

Mete:  Theorem  1 provides  a second  way  to  generate 
the  unit  homogeneous  Poisson  process  required  by  the 
time-scale  transformation  algorithm  previously  described. 
First  we  may  generate  a variate  from  the  Poisson  random 
variable  Nt  with  parameter  , say  Nt^  = n.  Then 
conditional  upon  = n,  n uniform  order  statistics  can  be 
generated  cn  the  interval.  For  reasons  explained  in 
Appendix  D,  this  second  method  was  used  in  the  computer 
program  implementation  of  the  time-scale  transformation 
method. 
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F.  RATIONALE  FOR  SELECTION  OF  DEGREE-TWO  EXPONENTIAL 
POLYNOHIAL  INTENSITT  FUNCTION 

The  intensity  function  identifying  the  non-honogeneous 
Poisson  process  investigated  in  this  paper  is  of  the  fora, 

X(t)  = exp(aQ  + a1t  + ct2t2) 

where  a^,  and  a 2 are  real  constants. 

This  intensity  function  was  selected  for  three  reasons. 
First,  an  intensity  function  aust  always  be  positive  (or 
zero)  if  it  is  to  be  meaningful.  The  above  function  is 
non-negative  for  all  values  of  and  Secondly, 

this  intensity  function,  by  proper  choice  of  constants,  can 
be  used  tc  represent  the  four  different  types  of  event 
streams  mentioned  previously  in  Section  I.  And  finally,  the 
selection  of  this  intensity  function  leads  to  simple 
statistical  procedures;  (for  details,  see  LEWIS  Ref.  9, 
p. 30-34) . 
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III. 


COMPETING  ALGORITHMS 


A.  GENEBAL 


Assuming  an  intensity  function  of  the  fora 


2 

X(t)  = exp(aQ  + + at  ) 

three  algorithms  for  generating  the  corresponding 
ncn-homogeneous  Poisson  process  are  discussed.  These 
algorithms  are  based  on  the  two  general  methods  presented  in 
Section  II-E,  and  the  decomposition  and  superposition 
property  of  Poisson  processes. 

B.  TIME-SCALE  TRANSFORMATION  ALGORITHM  (ALGORITHM  A) 


1 . Step  One 


Ey  definition  of  a non-homogeneous  Poisson  process, 
the  total  number  of  events  observed  over  a fixed  interval 
(0 , t«)  is  itself  a Poisson  distributed  random  variable,  N*.  , 

“ t0 

with  parameter  Ug  = Jg  X(t) dt.  The  first  step  cf  tne 
algorithm  is  to  determine  the  value  of  the  parameter  Ug  . 
Although  an  explicit,  closed-form  expression  for  the  above 
integral  cannot  be  found,  a series  representation  does 
exist.  Except  for  a constant  factor,  this  series 
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representaticn  assumes  the  fora  cf  the  error  function  cr  of 
Dawson's  integral.  A negative  value  for  the  coefficient  of 
the  second  degree  term,  a , yields  the  error  function  fora: 


whereas  a positive  value  for  results  in  the  Dawson's 
integral  fora: 


In  the  above  expressions  K,  t^,  and  t2  are  uniquely 
determined  by  the  coefficients  aQ,  a ct2  and  the  end  points 
of  the  interval  over  which  the  intensity  function  applies. 
A detailed  derivation  of  above  relationships  is  given  in 
Appendix  C. 

Evaluation  of  the  error  function  is  a FORTRAN 
supplied  procedure  and  requires  only  that  the  proper 
arguments  be  calculated  and  provided  to  the  FORTRAN  FUNCTION 
ERF  or  EERF  [Ref.  15].  Evaluation  of  Dawson's  integral  is 
best  accomplished  through  use  of  the  IMSL  (International 
Mathematical  and  Statistical  Libraries,  Inc.)  FUNCTION  MMDAH 
[Ref.  6].  The  accuracy  of  the  function  values  calculated  by 
these  routines  is  limited  only  by  the  precision 
characteristics  of  the  computer. 

2 . Ste£  Two 

Cnee  the  parameter  of  the  Poisson  random  variable 

N is  determined,  a realization  on  that  random  variable  is 
t0 
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required.  (The  approach  is  somewhat  backwards  since  it 
first  determines  hew  many  events  occurred  over  the  interval. 
It  then  distributes  that  fired  number  of  events  over  the 
interval  in  accordance  with  the  non-homogeneous  Pcisson 
process  described  by  the  intensity  function.  The  importance 
of  Theorem  1 now  becomes  evident  since  it  assures  the 
validity  of  such  a procedure.)  Generation  of  Poisson 

variates,  especially  those  with  large  parameter  values,  is  a 
complex  procedure  in  itself  if  efficiency  in  terms  of 
computer  time  and  memory  requirements  is  desired.  This 
problem  is  discussed  later  in  Section  IV-B.  For  the 
present,  assume  that  the  requisite  variate  has  been 
produced,  i.e.  H = n . 

fc0 

3 . Step  Three 

Given  that  n events  have  occurred  over  the  interval 
(0,tQ]  we  then  distribute  n events  along  an  interval  of 
length  UQ  in  accordance  with  a homogeneous  Poisson  process. 
Since  events  in  a homogeneous  Pcisson  process  are  uniformly 
distributed  over  an  interval  (given  that  n events  have 
occurred),  this  step  merely  requires  that  n uniform  (0,1) 
variates  be  generated,  ordered. from  lowest  to  highest,  and 
then  each  multiplied  by  the  factor  pg.  The  values  in  this 
n-element  vector,  (u^  , , . ..,  u^)  , correspond  to  the 

points  plotted  on  the  vertical  axis  in  Figure  3. 

4 . Step  Four 

Each  event  in  the  homogeneous  Poisson  process  must 
be  transformed  by  the  inverse  of  the  integral  of  the 
intensity  function.  Letting  ft  \ (s)  ds  = A (t)  , the  inverse 

i 0 

a “ (•)  applied  to  each  event  in  the  homogeneous  Pcisson 
process,  will  produce  a corresponding  event  in  the 
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non-homogeneous  Poisson  process,  i.e.  A = t.  , 

A”^(u2)  = t2#  etc.  The  difficulty  is  that  since  the 

integral  of  this  specific  intensity  function  cannot  usually 

be  explicitly  expressed,  the  fora  of  its  inverse  usually 

eludes  any  convenient  computational  formula  expression.  The 

unique  position  on  (0,tQ]  the  inverse  determines  for  each 

input  value  can  be  found  to  any  degree  of  accuracy  desired, 

by  iterative,  numerical  methods.  The  Newton-Raphson  method 

is  easily  employed  and  very  efficient  in  the  present 

scenario.  Its  implementation  is  explained  in  Section  IV-C. 

Since  the  function  A (t)  is  strictly  monotone  increasing, 

the  inverse  function  A~^  (u)  applied  to  an  ordered  sequence 

of  input  values  results  in  an  ordered  sequence  of  output 

values.  Therefore,  t,  , t , ...,  t are  the  times  cf  events 

12  n 

in  the  non-hcaogeneo us  Poisson  process  and  the  algorithm  is 
complete . 


C.  TIME-SCALE  TEA NSEORMATI ON  ALGORITHM,  ALT2BNATE 
(ALGORITHM  A') 


An  alternative  approach  to  the  time-scale  transformation 
method  described  above  is  to  generate  the  required 
homogeneous  Poisson  process  by  using  the  fact  that  in  this 
process  the  random  times  between  events  are  independently 
exponentially  distributed.  Thus  one  generates  unit 
exponential  variates  until  their  sum  exceeds  . The 
partial  sums  give  the  times  to  events  and  the  number  of 
partial  sums  less  than  or  equal  to  pQ  is  a Poisson  (Uq) 
variate.  Rote  that  the  Poisson  variate  comes  out  as  a 
by-product  in  this  procedure  rather  than  as  a pre-product  as 
in  Step  Two  of  Algorithm  A above. 

Although  this  method  combines  Step  Two  and  Step  Three  of 
Algorithm  A into  a single  procedure,  it  is  not  necessarily 
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the  best  method.  Because  it  requires  the  use  of  the 
additive  method  of  Poisson  variate  generation,  it  becomes 
inefficient  for  Poisson  processes  with  many  events  (see 
Appendix  D) . Por  this  reason  Algorithm  A instead  of 
Algorithm  A’  was  used  when  implementing  the  time-scale 
transformation  method  into  a computer  program. 

D.  POIS  SCN-DECOMPOSITION  AMD  GAP  STATISTIC  ALGOBIT  8M 
(ALGORITHM  B) 

s 

It  is  recommended  that  the  reader  refer  to  Pigures  4,  5 
and  6 for  a better  understanding  of  the  following  steps.) 

1 . Step  One 

The  Poisson-decomposit ion  and  gap  statistic 
algorithm  begins  with  an  examination  of  the  coefficients  of 
the  intensity  function.  By  doing  so,  the  intensity  function 
is  categorized  into  one  of  six  possible  configurations. 
These  six  cases  are  discussed  in  LEWIS  and  SHEDL2R 
[Ref.  11].  Examples  of  each  case  are  illustrated  in  Figures 
7 through  12  located  at  the  end  of  this  section. 


2 . Step  Two 

a.  If  the  intensity  function  \ (t)  is  monotone 
increasing  or  monotone  descreasing  over  the  interval  (Cases 
I,  II,  IV  and  V;  see  Figures  7,  8,  10  and  11)  the  intensity 
function  is  decomposed  into  two  separate  intensity  functions 
over  the  same  interval.  The  resulting  intensity  functions 
are  of  the  form; 

\(t)  = exp  (Yq  + Y1t) 

and 
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X*(t)  = X (t ) - A_(t)  = exp(aQ  + + a2t2)  - exp^  + Y^t) 

It  is  clear  that  X_ ( t ) + X*(t)  = X(t). 

b.  If  either  of  the  two  cases  not  covered  by  2a 
abcve  occurs,  X ( t)  will  be  monotone  increasing  (decreasing) 
on  the  subinterval  (0,b]  and  monotone  decreasing 
(increasing)  on  the  complementary  subinterval  (b,t0]  (see 
Figures  2 and  6)  , where  b is  a unique  point  within  the 
interval  at  which  X (t)  has  a maximum  (minimum)  value.  By 
dividing  the  interval  properly  into  two  disjoint,  contiguous 
subintervals,  each  subinterval  may  be  treated  as  explained 
in  2a.  Subsequent  steps  are  applied  to  each  of  the  two 
intervals  separately  and  the  results  combined. 

3 . Step  Three 

An  efficient  algorithm  for  generating  a 

ncn-homogeneous  Poisscn  process  with  a log-linear  intensity 
function  (i.e.  X ( t)  = exp  (3Q  ♦ fi^t)  ) is  presented  by  LEWIS 
and  SHEELEE  [Hef.  13].  This  algorithm  generates  the 

non-homogenecus  Poisson  process  through  the  use  of  gap 
statistics.  (A  comparison  of  the  gap  statistic  technique 
with  the  conventional  integral  transform  technique  is 

discussed  in  Appendix  C.)  By  judicious  selection  cf  the 
coefficients  cf  the  log- linear  intensity  function,  most  of 
the  total  area  under  the  original  intensity  function  x (t) 
will  be  contained  under  the  function  X (t) . Therefore  most 
of  the  events  in  the  non- homogeneous  Poisson  process  with 
intensity  function  X (t)  can  be  accounted  for  by  employing 
the  gap  statistics  algorithm  on  the  intensity  function 
X (t)  . 
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Step  3 


Generate  events  T.  from  log-linear  intensity 
function,  A_(t),  using  gap  statistic  algorithm. 
Events  will  be  ordered. 


1 2 
Generate  events  s^  from  intensity  function 
A*(t)  = x(t)  - Ajt)  using  rejection  algorithm. 
Events  will  not  be  ordered. 
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Step  4 
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Step  5a 

- Order  events  from  Step  4, 

T1  T2 

T3T4  T5  T6  T7  T8  T9 

T10 

, 1 1 

11  111  II 
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1 > 
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Step  5b 

- Merge  (superpose)  events  from  Steps 

3 and  ! 

Result  is  event  stream  T-j , T^,  ... 

\(t)  = \(t)  + A*(  t) . 

from 

Figure  5 - DIAGRAM  (CONTINUED) 
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4 . Ste p Four 

All  that  remains  to  be  done  is  to  generate  an 
ordered  sample  of  some  given  size,  on  the  interval  (0,tQ], 
from  the  remaining  component  of  the  original  intensity 
function,  i.e.  \*  (t)  . 

The  number  of  events  in  the  interval  is  a Poisson 

random  variable  Ni.  with  parameter  u'  = /t0x*(t)dt. 

0 

Once  a realization  on  this  random  variable  is  obtained,  i.e. 

Nf  = n'.  Theorem  1 may  be  invoked.  Since  A*(t)/y'  is 
t0  0 

more  easily  evaluated  than  its  indefinite  integral,  the 
rejection  technique  (explained  in  Section  IV-D)  is  used  to 
generate  the  n'  required  variates.  The  rejection  technique 
is  not,  in  general,  always  an  efficient  method  for  variate 
generation  unless  great  care  is  taken.  let  in  this 
deccmpositicn  scenario,  it  will  be  used  to  generate  only  a 
small  percent  of  the  total  events  required  by  the  original 
intensity  function  X(t) . The  majority  of  the  events  will 
be  generated  by  the  efficient  gap  statistics  algorithm.  The 
efficiency  gains  should  more  than  compensate  for  any 
efficiency  losses  due  to  the  use  of  the  rejection  technique. 

5 . Step  Five 


The  events  produced  by  Step  3 will  be  in  order  on 
the  interval  (0,t0],  The  events  produced  in  Step  4 will  not 
be  in  order  cn  (0 , tQ  ] . By  ordering  the  events  from  Step  4, 
(which  are  few  in  number  compared  to  the  total  number  of 
events  in  the  non-homogeneous  Poisson  process)  it  is 
possible  tc  superpose  the  two  ordered  event  streams.  The 
merged  event  streams  produce  a new  event  stream  from  the 
original  intensity  function  \ (t) . 
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Figure  11  - SAMPLE  INTENSITY  FUNCTION  - CASE  V 
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IV.  ALGORITHM  IMPLEMENTATION 


A.  GENERAL 

This  section  explains  the  several  specific  techniques 
that  were  applied  to  implement  the  two  competing  algorithms 
(Algorithm  A and  Algorithm  S)  into  FORTRAN  computer 
pregrams.  Detailed  discussion  of  the  various  subprograms  is 
avoided  since  References  11  and  13  and  attached  program 
listings  provide  such  information.  The  Algorithms  A and  B 
have  some  subprogram  requirements  in  common  while  other 
subprograms  are  unique  to  one  algorithm  or  the  other.  This 
section  will  discuss  first  those  requirements  common  to  both 
algorithms,  then  those  needed  only  by  the  time-scale 
transformation  algorithm  (Algorithm  A)  and  finally  those 
unique  to  the  Poisson-decomposition  and  gap  statistic 
algorithm  (Algorithm  B) . Hereafter  differentiation  between 
the  algorithms  and  the  FORTRAN  computer  program 
implementation  of  the  algorithms  will  not  always  be  made. 
The  meaning  will  be  clear  from  the  context. 

B.  COMMON  REQUIREMENTS 


1 . Integration  of  a Degree-Two  Exponential  Polynomial 
Function  over  a Fixed  Interval 

Algorithm  A requires  that  the  intensity  function 


l. 
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X (t)  be  integrated  over  a fixed  interval  in  order  to 
determine  the  value  of  the  parameter  of  the  Poisson  random 
variable  that  governs  the  number  of  events  that  are  observed 
to  occur  in  the  non-homogeneous  Poisson  process.  Algorithm  8 
requires  that  the  intensity  function  X*(t)  = A(t)  - X (t)  , 
be  integrated  over  a fixed  interval  for  the  same  reason. 
Since 

b b b 

/ A*<t)dt  = / X (t)  dt  - / X(t)dt 
a a a 

and  \ ( t ) has  an  explicit  expression  for  its  indefinite 
integral,  i.e. 

b 

/ X_(t)dt  = exp(Y0) 

cL 

the  problem  for  both  algorithms  is  reduced  to  computing  the 
value  for 

b b 2 

/ X(t)dt  = / exp(aQ  + a^t  + c^t  ) dt  . 
a a 

SUBROUTINE  HELP  employs  IMSL  FUNCTION  NMD AW  or  the  FORTRAN 
supplied  procedure  DERF,  as  appropriate,  to  perform  this 
calculation.  Section  III-3  and  Appendix  B discuss  how  this 
computation  could  also  be  made  using  a convergent  series 
representation. 


exp  (y^t ) 


2*  Generation  of  a Poisson  Varia  te  with  a Given 
Parameter 

Both  candidate  algorithms  require  at  least  one 
realization  on  a Poisson  distributed  random  variable  with  a 
given  parameter.  The  value  assumed  by  the  random  variable 
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is  the  number  of  events  observed  in  a specific  interval  of 
time  over  which  the  occurrence  of  events  is  governed  by  a 
given  intensity  function  (i.e.  X (t)  for  Algorithm  A and 
X * (t)  for  Algorithm  B)  . The  nature  of  most  real  event 
streams  which  lend  themselves  to  analysis  using  a 
non-homogeneous  Poisson  process  model  is  that  they  consist 
of  a large  total  number  of  events.  This  will  be  the  case 
either  if  a dense  process  is  observed  over  an  interval  of 
short  or  moderate  length  or  if  a '•sparse"  process  is 
observed  over  an  extended  interval  (see  the  arrivals  at  an 
intensive  care  unit  given  in  LEWIS  [Ref.  9]).  The  point  is 
that  both  algorithms  must  be  flexible  enough  to  generate 
ncn-homogeneous  Poisson  processes  that  result  in  high 
numbers  of  events  occurring.  Since  the  number  of  events 
occurring  over  any  interval  in  such  a process  is  a Poisson 
random  variable,  it  becomes  necessary  to  be  able  to  generate 
a variate  from  a Poisson  distribution  with  a large 
parameter,  i.e.  large  mean.  The  two  most  theoretically 
straightforward  algorithms  for  generating  Poisson 
distributed  variates  (the  additive  and  multiplicative 
methods)  become  computationally  cumbersome  as  the  parameter 
of  the  random  variable  increases.  Both  the  additive  and 
multiplicative  algorithms  and  the  practical  deficiencies  of 
each  are  discussed  in  Appendix  D. 

The  technique  selected  to  deliver  a Poisson  variate 
on  demand,  to  both  Algorithm  A and  Algorithem  B is  the  Gamma 
method.  It  is  developed  and  explained  in  an  unpublished 
book  by  AHRENS  and  DIETER  [Ref.  1],  A paraphrased  account 
of  the  development  of  this  algorithm  is  included  in 
Appendix  E.  The  main  advantage  of  AHRENS  and  DIETER'S  Gamma 
method  is  that  whereas  the  computer  time  required  for  the 
additive  and  multiplicative  methods  is  proportional  to  the 
parameter  of  the  Poisson  distribution  being  sampled,  the 
computer  time  required  by  the  Gamma  method  is  proportional 
to  the  logarithm  of  the  parameter.  The  SUBROUTINE  PVAR 
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employs  the  Gamma  algorithm  to  return  a Poisson  variate  when 
given  a parameter  as  an  input. 

3.  Hve^t  Storage 

Any  algorithm  which  simulates  a non-homogeneous 
Poisson  process  must  have  some  mechanism  to  provide  the  user 
with  information  that  completely  describes  the 
realization (s)  of  the  process  simulated.  Since  the  specific 
arrangement  of  the  stream  of  events  over  the  interval  is  the 
information  cf  interest  to  the  analyst,  the  location,  or 
time  of  occurrence,  of  each  single  event  on  the  interval 
must  be  stored.  Equivalently,  the  spacings  between  events 
would  completely  describe  the  realizations  on  a 
non-homogenecus  Poisson  process.  Thus  event  spacing 
information  rather  than  event  location  information  could  be 
stored.  {Programs  written  for  candidate  algorithms  A and  B 
both  provide  the  option  for  the  user  tc  demand  either  event 
location  or  event  spacing  information.)  when  using  either 
algorithm,  an  array  large  enough  to  hold  location  or  spacing 
data  for  each  event  generated  by  the  algorithm  must  be 
created.  Since  the  number  of  events  observed  in  any 
realization  of  a non-homogeneous  Poisson  process  is  itself  a 
random  variable,  the  precise  size  of  the  array  cannot  be 
determined  a priori . If  the  programs  are  to  have  any  value 
for  general  application,  they  must  be  able  to  accept 
intensity  function  parameters  which  will  demand  large 
numbers  of  events  when  simulated.  Using  the  somewhat 
arbitrary  assumption  that  an  event  stream  with  an  average  of 
4500  events  is  sufficient  for  most  simulation  scenarios,  a 
fixed  array  with  a capacity  for  5000  events  is  used  in  the 
pregrams  implementing  both  of  the  algorithms.  If  the  value 
of  the  intensity  function  integrated  over  the  interval  is  as 
high  as  the  maximum  limit  of  4500  (i.e.  the  number  of  events 
in  the  interval  is  Poisson  distributed  with  mean  = 4500) 
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then  the  array  of  length  5000  allows  the  Poisson  random 
variable  to  exceed  its  mean  by  7.45  standard  deviations 
before  an  array  overflow  is  encountered.  This  highly 
unlikely  event  is  of  such  rare  occurrence  (less  than  one 
chance  in  a billion)  that  it  may  be  disregarded.  However, 
should  it  cccur,  the  programs  will  reinitialize  themselves 
and  generate  a new  Poisson  process.  Also  an  error  indicator 
is  incremented.  This  error  indicator  (IER)  may  be  written 
on  demand.  Its  value  is  the  number  of  times  the  program  was 
forced  tc  abort  generating  a Poisson  process,  reinitialize 
itself  and  start  again. 

A 5000  element  array  is  small  enough  to  avoid  its 
being  an  undue  memory  requirement  burden  on  most  operating 
systems,  yet  large  enough  to  accommodate  most 
non-homogenecus  Poisson  processes  of  interest.  Its  choice, 
though  arbitrary,  was  based  upon  the  above  two 
considerations. 

C.  SPECIAL  REQOIREMENTS  OF  IHE  TIME-SCALE  TBANSPORMATICN 
ALGCBITHK  (ALGORITHM  A) 


1 . Onif era  Variates 

As  explained  in  Section  III-B,  once  the  number  of 

events  observed  in  the  non-homogeneo us  Poisson  process  is 

established  (i.e.  N.  = n)  , it  is  necessary  to  construct  a 

fc0 

homogeneous  Poisson  process  consisting  of  n events  over  an 
interval  of  length  units.  This  is  easily  done  by 

generating  n uniformly  distributed  variates  on  the  interval 
(0,1)  and  then  scaling  each  by  the  factor  . The  LLRANDOM 
computer  library  package  developed  by  LEWIS  and  LEABMONTH 
[Ref.  12]  is  a very  efficient  source  of  such  variates.  Once 
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an  array  of  n uniform  (0,1)  variates  is  obtained  from 
LLRANDOM , each  element  of  the  array  must  be  multiplied  by 
the  appropriate  scaling  factor.  It  is  these  resulting 
numbers  which  must  be  acted  upon  by  the  inverse  of  the 
integrated  intensity  function  to  yield  the  location  of 
events  in  the  non-homogeneous  Poisson  process  on  the 
original  interval  (0,tg]. 


2 • Sorting  of  Events 

Tc  he  of  much  practical  value  a simulation  routine 
for  a non-hcmogeneous  Poisson  process  should  order  the 
events  in  the  interval  from  first  to  last.  In  Algorithm  A, 
this  ordering  may  be  done  before  applying  the  inverse 
integrated  intensity  function  to  the  elements  in  the 
homogeneous  Poisson  process.  The  monotonic  nature  of  the 
integrated  intensity  function  maintains  the  relative  order 
of  all  elements  after  they  have  been  transformed  by  the 
inverse  function  (see  Figure  3) . Of  course  the  ordering 
could  be  dene  after  the  transformations  also.  In 
implementing  Algorithm  A,  the  uniform  variates  were  scaled 
by  the  factor  then  ordered,  from  lowest  to  highest, 
before  being  transformed  by  the  inverse  of  the  integrated 
intensity  function. 

Ordering  of  large  arrays  of  numbers  is  a time 
consuming  operation  on  the  computer.  There  are  many 
ordering  algorithms  of  varying  degrees  of  sophistication  and 
efficiency.  Because  ordering  is  unavoidable  when  using 
Algorithm  A,  selection  of  an  efficient  ordering  routine  is 
important.  The  ordering  algorithm  used  in  this 
implementation  was  the  W.  R . CHDRCH  computer  center  library 
routine  FXSOFT  which  employs  Singleton's  version  of  the 
partition  exchange  sort.  (A  program  listing  of  PXSCRT  is 
provided  in  the  computer  center's  catalogue  of  library 


51 


1 


routines.)  The  PXSORT  routine  appears  to  be  the  most 
efficient  ordering  routine  readily  available  for  the 
purpose.  (It  is  acknowledged  that  more  efficient  routines 
specifically  tailored  to  the  problem  of  ordering  uniform 
random  variates  may  possibly  have  been  developed  that  would 
improve  the  overall  efficiency  of  the  program  implementing 
Algorithm  A.) 

3.  Computation  of  the  Trans  formed  Values 

The  essence  of  Algorithm  A is  the  application  of  the 
inverse  cf  the  integrated  intensity  function  to  each  event 
in  the  homogeneous  Poisson  process  on  the  interval  (0,pq]. 
As  mentioned  before,  the  intensity  function,  A(t) , under 
consideration  does  not  yield  an  explicit  expression  for  the 
integrated  intensity  function,  A(t)  = ?(t),  that  must  be 
inverted.  (Note:  F(t)  is  a "scaled"  distribution  function). 
Hence  neither  does  a ccmputa tionally  convenient  expression 
for  this  inverse  function  exist.  Numerical  methods  must  be 
employed  to  transform  the  position  of  each  event  on  the 
interval  (0,Ug]  in  the  homogeneous  process  tc  its 
corresponding  position  on  the  interval  (0,tg]  over  which  the 
simulated  non-homogeneous  Poisson  process  is  to  be  produced. 

The  Newton-Ra phson  technique  was  used  to  accomplish 

the  transformation.  This  technique  allows  for  iterative 

approximations  which  converge  to  the  true  transformed 

values.  (i.e.  t^,  t2#  ...,  tn;  where  t ^ = F 1(ui))  The 

iterations  continue  for  each  value  u^_  until  an  estimate  t* 

for  its  corresponding  t^  is  found  that  satisfies 

IP(t’)  - u,  | < e , where  e is  a predetermined  tolerance 

-5 

level.  (The  specific  value  for  e used  was  uQ  x 10  .) 

In  the  Newton-Raphson  technique,  given  a function 
h (x) , the  objective  is  to  find  the  solution,  x*,  satisfying 
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the  equation  h (x*)  = 0.  Letting  x ^ be  the  initial  estimate 
for  x*  (based  upon  some  reasonable  criteria)  the  next 
estimate  for  x*,  that  is,  x2»  can  be  calculated  from  the 
fundamental  iteration  relationship, 

h(xk) 

Xk+1  = xk  h' (xk)  • 

This  assumes  that  the  function  h (•)  and  its  derivative  h * ( • ) 
can  be  evaluated  at  xk.  The  process  is  repeated  k times 
until  | h (xk) | <e  ; or,  until  | xk  - xk+1 | <e  . This 
procedure  is  nothing  more  than  using  the  first  two  terms  of 
the  Taylor  series  expansion  of  the  function  h(«)  to  locate 
the  x-intercept  of  the  line  tangent  to  h (x)  at  the  point 
h(xk).  That  intercept  value  becomes  the  next  approximation 
for  the  root  of  h (x)  and  the  procedure  is  repeated.  A 

graphical  explanation  of  the  Newton-Ra phson  method  is 
presented  in  Appendix  E. 

Applying  the  Newton-Raphson  method  to  the  present 

problem  requires  some  special  modifications.  The  problem  is 

to  find  a value  t.  such  that  F(t.)  = u. , where  u.  is  known, 

i 11  i 

i.e.  to  find  t.  = F (u. ) . Direct  application  of  the 

inverse  F ^(*)  is  impossible  since  its  functional  form 
defies  all  but  the  most  abstract  mathematical  expression. 
However  if  a new  function  G^  (t)  = F(t)  - u^  is  defined  and 
if  its  root  t*,  determined  (i.e.  a t*  such  that  G^  (t*)  - 0) 
th  • F(t*)  - u^  = 0,  and  F(t*)  = u^ . Thus  t^  = t*  is  the 
value  desired. 

In  order  to  apply  the  Newton-Raphson  method  to  the 
function  G ^(t)  it  must  be  possible  to  calculate  both  G^  (t) 
and  its  derivative.  Since  G^  (t)  differs  from  F (t)  only  by  a 
constant,  the  problem  reduces  to  that  of  evaluating  F (t) . 
This  has  already  been  done,  as  previously  explained  in 

Section  III-E,  and  merely  requires  the  assistance  of 
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SOEROOTINE  BELP.  Now  G^<t)  = F*  (t)  , so  taking  the 
derivative  cf  G^t)  returns  the  intensity  function  X(t) 
which  is  easily  evaluated  for  any  t.  Clearly,  the 
Newton-Raphson  method  nay  be  used. 

It  should  be  noted  here  that  for  each  u^  there  is  a 

corresponding  function  G^t)  on  which  the  Newton-Raphson 

technigue  must  be  used.  Since  each  use  of  the 

' Newton-Raphson  method  nay  require  several  iterations  to 

arrive  at  a value  t*  which  satisfies  the  tolerance 

criterion,  it  is  readily  evident  that  for  large  n, 

considerable  computational  effort  is  required  to  obtain  all 

the  values  t (i  * 1 , ...,  n)  . The  number  of  iterations 

required  to  arrive  at  a suitable  approximation  for  each  t. 

is  highly  sensitive  to  the  accuracy  of  the  initial 

approximation  (designated  t.  J If  the  initial  approximation 

I 

is  close  to  the  actual  t^,  the  Newton-Raphson  method  will 
converge  very  quickly.  If  the  initial  approximation  is 
poor,  convergence  could  be  much  slower.  The  procedure  used 
to  select  initial  approximations  for  each  t^  will  therefore 
have  a profound  effect  upon  the  overall  efficiency  of 
Algorithm  A. 

Selection  cf  these  initial  approximations  is  dene  by 
partitioning  the  interval  (0,tQ]  into  equal  length  segments. 
The  number  of  segments  is  equal  to  min[10,  n/4  ].  The 
function  F ( • ) is  evaluated  at  the  end  points  of  each  segment 
and  these  end  points,  with  their  corresponding  function 
values  are  stored  in  an  array.  This  procedure  is  performed 
by  SOBROOTIN2  BNCHMK . See  Appendix  E for  a graphical 
representation. 

To  find  an  initial  approximation  for  the  t^  which 
corresponds  to  any  given  u^,  the  array  of  function  values  is 
searched  until  two  adjacent  function  values  which  bracket 
are  located.  These  adjacent  function  values  uniquely 
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identify  the  segment  of  the  interval  (0,t^]  in  which  the 
elusive  t^  is  located.  Either  end  point  cf  the  segment 
would  serve  as  a good  initial  approximation  t^^fcr  the 
Newton-Haphson  method.  However,  for  the  purpose  of  this 
implementation,  the  end  point  which  yielded  the  function 
value  closer  to  u^  is  used  for  the  initial  approximation. 

The  decision  to  divide  the  interval  (0,tQ]  into 
min  [10,  n/4  ] segments  was  based  upon  empirical  results.  Of 

t 

several  proposed  segmenting  schemes  that  cne  which  resulted 
in  the  fewest  total  number  of  distribution  function 
evaluations  over  several  intensity  function/interval 
scenarios  was  chosen.  Higher  values  than  n/4  for  sparse 
event  streams  may  produce  better  results.  But  n/4  gave  good 
results  for  dense  streams  and  adequate  results  for  sparse 
streams.  It  appeared  to  be  the  best  compromise  as  a 
candidate  for  general  usage.  (The  number  of  function 
evaluations  of  G^(t)  for  each  t^  averaged  between  2.2  and 
2.7  for  this  partitioning  scheme.) 

Alternative  methods  would  be  to  use  for  the  initial 
approximation  for  t one  of  the  following  values: 

• *1-1  «*i> 

" *1-1  * t/n 

Neither  cf  these  other  methods  were  attempted  so  it  is 
uncertain  whether  they  would  be  more  or  less  efficient  than 
the  method  which  was  used.  Efficiency  differences  ataong  the 
three  methods  are  probably  not  substantial  since  each  will 
usually  yield  a good  first  approximation. 
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D.  SPECIAL  REQUIREMENTS  OP  THE  DECOMPOSITION  ALGORITHM 
(ALGORITHM  B) 


1 . Intensity  Function  Categorization 


The  decoaposition  routine  selected  depends  on  which 
of  six  possible  shape  categories  the  intensity  function 
A (t)  falls  into  (c.f.  Figures  7 through  12)  . An 
exaaination  cf  the  constants  aQ , and  a2  the  intensity 
function  and  the  interval  (0,tg]  over  which  the 
non-homogeneous  Poisson  process  is  to  occur  will  uniquely 
identify  the  category  of  the  intensity  function.  A thorough 
discussion  of  this  procedure  is  presented  in  Ref.  11,  and  is 
net  reproduced  here.  implementation  of  the  category 
identification  procedure  requires  a lengthy  sequence  of 
decision  statements  within  the  computer  program. 

2.  Selection  of  the  Imbedded  Log-Linear  Intensity 
Function  (s) 


Once  the  intensity  function  has  been  categorized  it 
must  be  decomposed  in  accordance  with  a complex  scheae 
described  in  Ref.  11.  The  objective  of  the  decomposition 
scheme  is  tc  fit  a leg-linear  curve  (or  curves)  coapletely 
underneath  X (t)  on  the  interval  (i.e.  X^(t)<  X (t)  , 0<t<tg) 

in  such  a way  as  to  maximize  the  area  under  the  log-linear 
curve  (s).  This  is  done  by  partitioning  of  the  interval 
(0,tQ]  into  subintervals  (0,b]  and  (b,tQ]  if  necessary,  and 
then  by  proper  selection  of  the  coefficients  yq  and  y^  for 
the  log-linear  f unction (s).  These  coefficients  are 

functions  cf  the  coefficients  aQ , and  in  the  original 
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intensity  function  X (t)  and  of  the  interval  (0,tg].  As  the 
program  advances  through  the  categorizing  decision 
statements,  the  proper  coefficients  are  computed  for  the 
appropriate  leg-linear  function(s)  . 


3.  Gag  Statistics  Algor ithm 


The  precursor  to  L2WIS  and  SHEDLER  [Hef.  11]  was  a 
paper  by  the  sane  authors  [Ref.  13]  proposing  a gap 
statistics  algorithm  for  simulating  a non-honogeneous 
Poisson  process  with  a log-linear  intensity  function 
\ (s)  = exp  (Sq  ♦ g^s) . (Reference  11  reviews  this  algorithm 
in  detail.)  As  previously  noted.  Algorithm  B divides  the 
intensity  function  \ (t)  into  a sum  of  two  intensity 
functions,  Mt)  and  X«(t).  The  new  intensity  function 
Mt)  is  chosen  to  take  on  the  fora  of  a log-linear  function 
so  that  the  cap  statistics  algorithm  may  be  used  to  generate 
a stream  of  events  from  this  portion  of  the  original 
intensity  function  X (t) . The  SUBROUTINE  NHPP2  implements 
the  LEWIS  and  SH2DLER  gap  statistics  algorithm  for  the  input 
intensity  function  X^(t)  and  returns  an  appropriate  event 
stream  to  the  calling  program.  Since  the  use  of  the  gap 
statistics  algorithm  within  Algorithm  3 is  the  rationale  for 
suspecting  it  to  be  a more  efficient  algorithm  than  the 
time-scale  transformation,  it  is  instructive  to  examine  the 
efficiencies  gained  in  using  the  gap  statistics  algorithm 
vice  a time-scale  transformation  algorithm  on  a leg-linear 
intensity  function.  This  question  is  addressed  in  Appendix 
C.  It  was  found  that  use  of  the  gap  statistics  algorithm 
resulted  in  a reduction  in  computer  time  of  approximately 
50*  from  that  required  by  the  time-scale  transformation 
method . 
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4 • The  Be  lection  Routine 

The  gap  statistics  algorithm  has  efficiently 
produced  an  event  stream  from  a non-hoaogeneous  Poisson 
process  defined  by  the  log-linear  intensity  function 
X (t)  = exp  (y0  ♦ y-^t)  . But  the  process  to  be  simulated  has 
intensity  function  x (*)  * exp  (aQ  ♦ a^t  ♦ a2t2).  It  is 
therefore  necessary  to  superpose  an  event  stream  frcm  the 
intensity  function 

XMt)  = X(t)  - X(t) 

2 

* exp(aQ  + a^t  + a2t  ) - exp(yQ  + y^t) 

onto  the  event  stream  obtained  from  the  log-linear  intensity 
function. 

Once  it  is  determined  how  many  events  are  to  occur 
over  an  interval  with  intensity  function  X*(t),  i.e. 

= n',  it  is  necessary  tc  select  an  ordered  sample  of  n' 
variates  from  a distribution  with  density  function 


/ u X*(t)dt 
0 

The  value  n'  will  be  small  compared  to  the  total  number  of 
events  in  the  original  non-homogeneous  Poisson  process. 
Therefore  the  need  for  realizing  efficiency  in  generating 
these  n*  ordered  variates  is  not  usually  as  crucial  to 
overall  algorithm  performance  as  is  the  need  for  efficiency 
in  producing  the  non-homogeneous  Poisson  process  from  the 
log-linear  intensity,  X (t) , in  Step  3.  However  if  the 
technique  for  obtaining  these  n'  ordered  variates  is 
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extremely  inefficient,  much  or  all  of  the  efficiency  gains 
realized  frcm  Step  3 above  cculd  be  lost  here.  (In  fact  an 
example  of  such  a loss  of  previously  gained  efficiency  is 
documented  in  this  thesis;  see  Section  VI-B.) 

The  rejection  method  is  particulary  useful  for 
generating  random  variates  from  populations  with  continuous 
densities  that  are  bounded  and  which  are  concentrated  on  a 
finite  interval;  this  is  the  case  for  f*(t) . The  rejection 
method  and  an  algorithm  for  its  implementation  are  presented 
in  LEWIS  and  SHEDLER  [Ref.  11].  A geometric  argument  will 
suffice  to  exhibit  the  principle  involved.  Consider  the 
density  function  f (x) , for  a random  variable  X on  the 
interval  (a,b)  in  Figure  13.  The  maximum  ordinate  of  this 
density  is  c.  If  another  function  g (x)  * c*  > c is 
constructed  then  the  density  function  f(x)  is  enclosed 
within  the  rectangle  (a,0)  , (b,0),  (a  ,c*)  , (b,c*)  . If  a 

point  within  this  rectangle  is  selected  at  random  it  will 
fall  either  under  the  density  function  or  above  it.  If  it 
is  under  the  density  curve  the  abscissa  of  that  point  is 
accepted  as  a variate  from  the  population.  If  the  point 
lies  above  the  density  curve  that  point  is  rejected  and 
another  point  within  the  rectangle  is  selected  at  random. 
The  procedure  continues  until  n*  variates  have  been 
produced.  The  random  points  within  the  rectangle  are  easily 
produced  by  generating  two  independent  uniform  (0,1) 
variates  and  scaling  them  properly  resulting  in  a point 
(x,y)  , where  x * a ♦ u^(b  - a)  and  y = U2#c*.  Then  if 
fU)  * y,  x is  accepted  as  a variate,  otherwise  it  is  not. 
After  n'  variates  have  been  obtained,  they  are  ordered  to 
yield  an  event  stream  for  a process  with  intensity  function 
X*(t).  SOEROtJTINE  REJECT  employs  this  method  in  the 
program  for  Algorithm  B. 
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The  density  f(x)  is  comDletely  enclosed  by  the  rectangle 
(a,0),  (b,0),  (a,c*) , (b,c*). 

Procedure: 

1.  Generate  2 independent  uniform  (0,11  variates  u^ 
and  u 

2.  Compute:  x = a + bij^ , y = 

3.  Plot  (x,y). 

4.  If  the  point  (x,y)  lies  under  the  density  f(x)  accept  x 
as  a variate;  otherwise  go  to  Step  1. 

Example : In  the  graphical  example  above  x^  would  be  accepted 
as  a variate  whereas  would  not  be  accepted. 
(Note:  Ideally,  c*  = c and  a and  b correspond 
to  the  lateral  limits  of  the  density  *(x).  This 
minimizes  '"ejection  region.) 


.gure  13  - THE  REJECT I ON -ACCEPTANCE  METHOD  OF  7A9IATE 
GENERATION  PROM  AN  AR3ITR  ARY  DENSITY 


Since  the  probability  that  a point  on  the  abscissa 
will  not  be  accepted  as  a variate  is  proportional  tc  the 
area  within  the  rectangle  that  is  not  under  the  density 
function,  it  is  advantageous  to  Bake  this  area  as  saall  as 
possible.  Therefore  it  is  best  to  set  c*  = c and  set  the 
values  a and  b equal  to  the  end  points  of  the  interval  over 
which  the  density  function  occurs.  This  is  desirable  but 
not  necessary.  Figure  13  illustrates  the  sore  general  case 
where  the  rectangle  is  larger  than  necessary.  One  may  be 
required  to  select  a c*  > c if  c is  difficult  to  determine. 
Seldcn  would  a and  b not  coincide  with  the  lateral  liaits  of 
the  density  however. 

It  is  obvious  from  Figure  13  and  the  preceding 
discussion  cf  the  rejection  method  of  variate  generation 
that  it  is  the  "shape"  of  the  density  function  which  is 
critical  to  the  validity  of  the  method  and  not  the  fact  that 
it  integrates  to  one.  Therefore  any  scaled  density  would 
preserve  the  relative  shape  of  the  density  function.  As 
long  as  the  value  c*  is  at  least  as  great  as  the  aaxiaua 
value  of  the  scaled  density,  the  rejection  method  will  yield 
valid  variates.  The  intensity  f unction  X * (t)  may  be  thought 

of  as  a density  scaled  by  the  factor  y'  = /tQX*(t)dt. 

0 0 

Algorithm  B uses  the  intensity  function  X*  (t)  rather  than 
the  density  function  in  the  rejection  routine. 
This  eliminates  a division  step  for  every  function 
evaluation  required  by  the  method. 

5.  jerqinq  Event  Streams 

The  event  streams  from  the  intensity  functions  \ (t) 
and  X*  (t)  must  be  merged  to  produce  an  event  stream  from 
X (t)  = x(t)  ♦ x * (t)  • In  the  present  case  there  are  two 
event  streams  (one  long  and  one  short)  which  are  already 
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ordered.  This  structure  allows  for  superposition  rather 
than  a more  general  (and  more  tine  consuming)  ordering 
procedure.  To  do  this  merging  job  SUBROUTINE  COLATE  is 
called. 


E.  SUMMARY 


This  section  has  presented  a general  discussion  of  the 
more  significant  and  interesting  procedures  used  to 
implement  Algorithms  A and  B efficiently  in  FORTRAN.  The 
reader  is  referred  to  LEWIS  and  SCHEDLER  [Ref.  11]  for  a 
detailed  listing  of  steps  for  Algorithm  B.  Program  listings 
may  also  he  consulted.  These  may  be  found  after  Appendix  5. 


V.  METHOD  OP  COMPARISON  OF  ALGORITHMS 


A.  GENERAL 

Conclusions  drawn  from  the  results  of  comparing  the 
efficiencies  of  two  computer  programs  are  only  as  sound  as 
the  measures  of  effectiveness  upon  which  they  are  based. 
These  measures  of  effectiveness  are  always  somewhat 
arbitrary,  so  their  selection  should  be  based  upon 
compelling  reasoning. 

Even  after  reasonable  methods  of  effectiveness  have  been 
defined,  only  a finite  number  of  trials  (usually  a small 
finite  number)  for  a given  set  of  circumstances  may  be  run. 
Therefore,  general  statements  such  as  ’’this  algorithm  is 
better  than  that  algorithm"  can  seldom  be  made  with  complete 
confidence.  Yet  if  the  incomplete  information  gathered 
reveals  certain  trends,  generalizations  hypothesized  from 
such  trends  may  warrant  a high  degree  of  confidence. 

This  section  describes  how  the  two  algorithms  are 
compared  and  the  rationale  for  choosing  the  method  employed. 
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B.  HEASUEES  OF  EFFECTIVENESS 


1 . • Computational  Speed 

Perhaps  the  most  popular  efficiency  measure  for 
computer  programs  is  their  speed  of  execution.  Speed  is  a 
natural  standard  as  long  as  computer  time  remains  a costly, 
scarce  resource. 

Algorithms  A and  B were  compared  in  terms  of  the 
mean  time  required  to  generate  event  streams  from  a given 
set  of  intensity  functions.  It  should  be  noted  that  since 
the  number  of  events  in  each  event  stream  is  a random 
variable,  the  time  required  to  generate  one  event  stream  is 
also  a random  variable. 

Several  event  streams  with  different  intensity 
functions,  each  having  a different  expected  number  of 
events,  are  produced.  Event  streams  for  each  of  these 
intensity  functions  are  replicated  many  times  by  each 
algorithm.  The  total  execution  time  required  for  all 
replications  is  the  "speed"  measure  of  effectiveness. 

The  number  of  replications  varies  with  each 
intensity  function.  For  example,  an  event  stream  with  an 
expected  number  of  events  less  than  200  was  replicated  100 
times  whereas  event  streams  with  over  1000  expected  events 
were  replicated  only  30  times.  In  both  cases  the  total 
number  of  events  produced  was  of  the  same  order  of 
magnitude.  A priori  it  seemed  reasonable  to  assume  that 
computer  time  required  to  execute  each  algorithm's  program 
would  be  roughly  proportional  to  the  total  number  of  events 
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generated.  Therefore  in  designing  the  experiment  the 
product  cf  the  expected  number  of  events  times  the  number  of 
replications  (E[Ntg]  x # reps)  was  kept  roughly  constant  in 
order  to  keep  demands  on  computer  time  reasonably  under 
control.  (Note:  Program  execution  times  were  isolated  from 
compilation  and  linkage  times  for  the  purpose  of  the 
computational  speed  comparison.) 


2.  Computer  Memory  Requirements 


A second  natural  efficiency  measure  for  computer 
programs  is  their  core  memory  requirement.  Both  programs 
were  written  with  the  goal  in  mind  of  conserving  as  much 
memory  as  possible  without  unduly  increasing  program 
complexity.  Cere  requirements  reductions  could  most  likely 
be  accomplished  in  both  programs  through  more  sophisticated 
programming  methods.  Therefore  this  comparison  has  a major 
weakness  as  a measure  of  effectiveness.  Recognizing  this 
caveat,  memory  requirements  for  each  program  were  compared. 


3 . Fidelity 


The  question  of  paramount  importance  is:  How  well 
does  the  program  simulate  the  true  process?  This  is  a 
difficult  question  to  answer  when  dealing  with  a finite 
number  cf  realizations  on  a process  whose  theoretical 
properties  can  be  validated  only  after  an  infinite  number  of 
realizations. 


Several  methods  of  comparison  are  possible.  The 
simplest  is  that  cf  comparing  the  sample  means  and  sample 
variances  of  the  random  variable  whose  realizations  are  the 
number  of  events  generated  on  an  interval  by  a given 
intensity  function.  This  random  variable  should  be  Poisson 
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distributed  with  mean  equal  to  the  integral  of  the  intensity 
function  over  the  interval.  For  a Poisson  random  variable, 
the  variance  equals  the  mean.  This  test,  although  useful 
for  determining  if  either  algorithm  varies  grossly  from 
expected  performance,  yields  no  information  about  hew  the 
algorithms  are  distributing  their  simulated  events  over  the 
interval.  It  is  necessary  to  look  at  discrete  segments  of 
the  interval  and  examine  the  distribution  of  the  number  of 
events  produced  in  each  segment  by  our  simulated 
non-homogenecus  Poisson  process. 

From  the  definition  of  a non-homogeneous  Poisson 
process  (see  Section  II-B)  we  know  that  the  number  of  events 
observed  over  any  segment  of  the  interval  must  be  a Poisson 
random  variable  with  parameter  equal  to  the  integral  of  the 
intensity  function  evaluated  over  the  interval.  3y  dividing 
into  several  segments  the  interval  over  which  the  Pcisson 
process  is  being  simulated,  the  number  of  events  in  each 
segment  can  be  observed.  Two  hypothesis  tests  may  then  be 
made  for  each  segment. 


The  first  test  is  the  dispersion  test  [Ref.  3, 

p.  158].  This  test  seeks  to  ascertain  if  the  number  of 

events  observed  over  each  segment  is  distributed  as  a 

Poisson  random  variable.  Let  n,  , n_,  . ..,  n,  be  k 

12  k 

observations  on  the  number  of  events  occurring  in  a given 
fixed  segment  p,  and  let  n^=  (n^+  ...+n^)/k.  If  is  a 
Poisson  distributed  random  variable,  then  the  statistic 
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has  a distribution  which  is  approximated  by  a chi-sguared 
distribution  with  k - 1 degrees  of  freedom.  If  the  interval 
has  been  partitioned  into  m segments,  then  by  simulating  the 
nen-homogeneous  Poisson  process  k times,  we  can  perform  m 
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hypothesis  tests  at  a selected  significance  level  a . The 

test  statistic  dp  would  then  be  compared  to  the  (1-a) 

percentile  of  the  chi-squared  distribution  with  k - 1 

degrees  cf  freedom.  For  a large  number  of  hypothesis  tests, 

say  1000,  we  would  expect  approximately  1000a  rejections  of 

the  null  hypothesis,  if  in  fact  the  N 's  are  Poisson 

P 

distributed. 

For  example  consider  the  non-homogeneous  Poisson 
process  over  the  interval  (0,100].  If  this  interval  were 
partitioned  into  100  unit  segments  the  number  of  events 
occurring  in  each  segment  should  be  Poisson  distributed.  If 
51  event  streams  are  simulated  over  the  interval,  a value 
for  the  test  statistic  d may  be  computed  for  each  segment, 
i.e. 


Assuming  a significance  level  a-  .05  and  comparing  each  d 

2 P 

to  the  critical  value  x50  a=  67.5048,  we  would  count  the 

number  cf  test  statistics  such  that  dp  > 67.5048.  if  the 

number  of  events  in  the  segments  are  indeed  Poisson  random 

variables  we  would  expect,  on  the  average,  that  5 out  cf  the 

100  dp  values  would  exceed  the  critical  value. 

In  the  above  discussion  of  the  dispersion  test,  no 
mention  was  made  of  the  parameters  of  the  Poisson  random 
variables  (assuming  that  they  are  Poisson) . This  is  because 
the  dispersion  test  for  Poisson  distributions  is  parameter 
independent.  Thus  it  is  necessary  to  perform  a second 
hypothesis  test  to  determine  if  the  mean  number  of  events 
per  segment  is  equal  to  the  difference  of  the  integrated 
intensity  function  evaluated  at  the  right  and  left  end 
points  of  the  segment  respectively.  For  large  k the  sample 
mean  n is  normally  distributed.  Under  the  null  hypothesis. 
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a. 

i.e.  that  the  mean  of  ND  = /Dx(t)dt  * Ur>  (where  a.  and  a.  are 

^ P 3-3 

the  end  feints  of  segment  p)  n is  normally  distributed  with 

mean  Up  and  variance  Up/k.  The  test  statistic  is 


and  the  critical  value  is  the  (1-2. ) percentile  of  the 
normal  distribution,  za/2*  Again,  there  would  be  one 
hypothesis  test  for  each  segment  and  if  we  made  1000  such 
tests  we  would  expect  to  reject  the  cull  hypothesis 
approximately  1000a  times,  on  the  average. 

Continuing  the  example  above,  it  is  now  of  interest 
to  see  if  the  mean  number  of  events  on  each  segment  agrees 
with  the  theoretical  value  determined  by  the  integrated 
intensity  function.  This  time  assuming  a significance  level 
of  a = .10,  the  critical  value  for  each  test  statistic  Cp 
would  be  z^2=1.6<5.  Each  cp  > 1.645  would  be  counted  and 
if  in  fact  tie  random  variables  have  the  hypothesized  means, 
we  would  expect  (on  the  average)  that  approximately  10  test 
statistics  cut  of  100  would  exceed  1.645. 

After  the  dispersion  test  and  hypothesis  test  for 
the  means  have  been  performed  on  event  streams  from  each 
algorithm,  the  test  results  may  be  compared  to  see  if  either 
algorithm  appears  to  have  null  hypothesis  rejection 
percentages  which  are  not  consistent  with  the  a levels  of 
the  hypothesis  tests. 

(It  should  be  noted  that  during  the  development  of 
the  programs  histograms  of  the  event  streams  were  plotted  to 
determine  whether  or  not  they  ''fit”  their  respective 
intensity  functions.  Both  programs  seemed  to  produce 
satisfactory  results  by  this  subjective  measure.) 
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C.  OTHER  CONSIDERATIONS 


1 • Intensity  Function  Category 

Algorithm  B classifies  the  input  intensity  function 
into  one  of  six  categories  and  then  determines  what  action 
to  take  tc  generate  the  appropriate  event  stream.  It  might 
be  expected  that  the  speed  of  Algorithm  B will  be  affected 
not  only  by  the  total  number  of  events  to  be  generated  but 
also  by  the  form  of  the  input  intensity  function. 
Decomposition  of  the  intensity  function  into  four  sections 
will  prctably  require  more  time  than  if  it  must  be  divided 
into  just  two  sections. 

Algorithm  A treats  all  intensity  functions  alike. 
Therefore,  it  is  important  to  compare  the  two  algorithms  for 
all  six  possible  categories  of  intensity  functions.  It  is 
possible  that  Algorithm  B could  be  faster  than  Algorithm  A 
for  one  category  of  intensity  function  and  slower  than 
Algorithm  A fcr  a different  category  of  intensity  function. 
The  six  intensity  functions  selected  for  the  comparison 
represent  all  of  the  categories  defined  by  Algorithm  B,  and 
are,  in  fact,  those  intensity  functions  illustrated  in 
Pigures  7 through  12  of  Section  III.  These  categories  are 
numbered  from  I through  VI  and  correspond  directly  to  those 
listed  in  Ref.  11. 

2.  Evaluation  of  c«  Value  in  Rejection  Routine 

The  importance  of  reducing  the  size  of  the  rectangle 
enclosing  the  intensity  function  X*  (t)  prior  to  beginning 
the  rejection  algorithm  was  mentioned  in  Section  IV-D.  If 
possible  c*  should  correspond  to  the  maximum  value  of  A*(t) 


69 


r 


on  the  interval  (0,tQ].  Por  Cases  I and  II  (Pigures  7 and  8) 
this  maximum  value  is  easily  determined  since  it  occurs  at 
one  of  the  end  points  of  the  interval.  For  Case  III  (Figure 
9)  the  maximum  values  of  x*(t)  and  x£(t)  both  occur  at  the 
point  b where  the  interval  (0,tQ]  has  been  partitioned  into 
two  subintervals  (0 , b ] and  (b,tQ  ].  Por  Cases  IV,  V and  VI, 
(Pigures  10,  11  and  12)  the  maximum  value  of  X*(t)  is  not  so 
easily  determined.  LEWIS  and  SHEDLER  [Ref.  11,  p.  11]  state 
that  it  is  not  possible  to  find  this  maximum  explicitly. 
However  an  upper  bound  may  be  determined  from  the  values  k 
and  k where  X(t)  £ k and  x (t)  > it  on  the  interval  (0,*q  ] (or 
subintervals  (0,b]  and  (b,tQ  ] for  Case  VI).  Then  by  setting 
c*  = k - k it  is  certain  that  X*(t)  < c*  on  the  interval  or 
sutintervals  cf  interest.  Figure  14  which  illustrates  the 
rejection  routine  for  the  intensity  functions  for  Case  V and 
VI  used  in  the  analysis,  reveals  that  the  c*  computed  in 
such  a manner  may  not  be  very  efficient  as  an  upper  bound 
for  the  rejection  algorithm.  One  would  expect  some 
reduction  in  speed  efficiency  for  Algorithm  B when  c* 
greatly  exceeds  the  maximum  value  of  X*(t). 

3.  Designated  Tolerance  Level 

Algorithm  A requires  selection  of  a tolerance  level 
which  determines  the  stopping  criterion  for  Newton-Raphson 
iterations.  The  speed  of  the  algorithm  could  be  expected  tc 
be  very  sensitive  to  this  tolerance  limit.  By  setting  it 
tco  small  the  comparison  would  be  prejudiced  in  favor  of 
Algorithm  B.  By  not  setting  it  small  enough  the  fidelity  of 
the  resulting  event  stream  to  the  true  intensity  function 
would  be  impaired. 

It  was  determined  that  a tolerance  level  of 
X occasionally  demanded  that  the  computer 

discriminate  beyond  its  significant  digit  capability  in 
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single  precision.  This  resulted  in  the  failure  of  the 

Newton-Raphson  iterations  to  converge,  although  in  theory 

they  would  have  done  so  if  enough  precision  had  been 

retained.  A tolerance  level  of  B[H.  ] x 10~6  seemed  to  be 

c0 

unnecessarily  demanding  on  Algorithm  A even  though  the 
computer  couxd  discriminate  at  this  level.  A stopping  rule 
fixing  £ = E[Nt^]  x 10  ^ was  the  compromise  that  insured 
reasonable  accuracy  of  the  simulated  event  stream  without 
unfairly  burdening  Algorithm  A. 
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VI.  RESULTS  AND  CONCLUSIONS 


A.  GENERAL 


Computer  programs  implementing  Algorithms  A and  E were 
compared  in  terms  of  efficiency  as  explained  in  Section  V. 
This  section  documents  the  results  of  the  comparison  and 
discusses  general  observations  of  the  author  concerning  the 
two  competing  programs.  Also  recommendations  for  further 
study  concerning  Algorithm  B are  offered. 


B.  MEASURES  Of  EFFECTIVENESS  RESULTS 


1 . Speed 


Table  I convincingly  illustrates  the  superiority  of 
Algorithm  B over  Algorithm  A when  computational  speed  is 
used  as  a measure  of  effectiveness.  The  column  in  Table  I 
labeled  "Time  A/Time  B"  is  the  ratio  of  the  total  time 
reguired  by  Algorithm  A to  produce  a specified  number  of 
replications  of  the  given  non-homogeneous  Poisson  process, 
to  the  time  required  by  Algorithm  3 to  produce  the  same 
number  cf  replications  of  the  process.  For  all  six 
representative  intensity  functions  the  P cisscn-decomposition 
and  gap  statistic  algorithm  is  at  least  twice  as  fast  as  the 
time-scale  transformation.  For  large  event  streams  the 
superiority  cf  Algorithm  B is  even  more  pronounced. 
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Table  I - COMPUTATION  TIMES  FOR  EVENT  STREAMS 
(IBM  Systen  360/67,  FORTRAN  IV  (G)  con-pil^r) 


Results  obtained  for  2nd  trial  with  alternate  c*  value  determined  by  graphical  analysis. 


The  most  obvious  disparity  revealed  in  Table  I is 
the  difference  in  speed  efficiency  (of  Algorithm  E with 
respect  to  Algorithm  A)  between  the  group  of  Cases  I,  II  and 
III;  and  the  group  of  Cases  IV,  Y and  VI.  In  the  former 
group  Algorithm  B was  5 to  7.6  times  as  fast  as  Algorithm  A. 
For  the  latter  group  Algorithm  B was  only  2.2  to  2.5  times 
as  fast  as  Algorithm  A. 

This  difference  between  the  two  groups  immediately 
made  suspect  the  value  c*  used  in  the  rejection  routine  of 
Algorithm  B.  Por  the  first  group  c*  is  set  at  the  least 
upper  bound  for  \* (t) . For  the  second  group  c*  is  not,  in 
general,  a least  upper  bound  for  A*(t)  but  is  only  an  upper 
bound.  To  determine  whether  or  not  this  difference  in  the 
"guality"  of  the  c*  values  could  have  such  an  adverse  effect 
on  time  efficiency  for  Cases  IV,  V and  VI,  these  cases  were 
simulated  a second  time  with  new  c*  values  which  were 
empirically  evaluated  by  graphical  methods.  These  modified 
c*  values  were  set  close  to  the  least  upper  bounds  of  the 
respective  A*(t)  component  of  the  intensity  functions  for 
each  case.  Marked  time  efficiency  gains  were  observed  (see 
Discussion  column.  Table  I)  indicating  that  significant 
efficiency  losses  can  be  sustained  due  to  the  method  of 
setting  c*  values  in  Cases  IV,  V and  VI. 

* 

Three  factors  appear  to  influence  the  degree  of 
speed  efficiency  gains  of  Algorithm  3 over  Algorithm  A. 
First  is  the  selection  of  the  c*  value  as  discussed  above. 
The  second  factor  is  the  percentage  of  total  events  which 
Algorithm  E must  generate  from  the  rejection  routine.  For 
the  specific  intensity  functions  analyzed,  this  percentage 
ranged  from  4%  in  Case  V to  24S  in  Case  VI.  The  fewer 
events  required  from  the  rejection  algorithm  relative  to  the 
number  generated  from  the  gap  statistics  algorithm,  the 
greater  the  efficiency  of  Algorithm  3. 
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Finally  the  total  number  of  events  in  an  event 
stream  affects  relative  efficiencies  between  Algorithms  A 
and  B.  This  is  because  as  the  number  of  events  to  be  sorted 
by  Algorithm  A grows,  the  less  efficient  its  sorting  routine 
becomes.  Sorting  requirements  for  Algorithm  B are  greatly 
reduced  due  to  the  properties  of  the  gap  statistic  procedure 
used.  Therefore  as  total  events  increase,  so  does  the 
efficiency  advantage  of  Algorithm  B over  Algorithm  A 
increase . 


2.  Computer  Memory  Requirements 

Algorithm  A requires  approximately  72,000  bytes  of 
core  storage.  Algorithm  3 demands  about  92,000  bytes.  Both 
programs  are  costly  in  terms  of  storage  due  to  the  number  of 
library  routines  used  and  the  need  to  store  up  to  5000  event 
times.  The  difference  of  20,000  bytes  would  not  be 
important  unless  computer  capacity  were  being  challenged, 
tlemory  requirements  for  both  algorithms  could  be  reduced  by 
reducing  the  event  stream  capacity  of  the  programs.  A 
reduction  of  capacity  from  a maximum  of  5000  events  to  a 
maximum  cf  2000  events  in  each  program  would  result  in 
storage  requirements  of  54,000  bytes  and  68,000  bytes  for 
Algorithm  A and  Algorithm  B respectively. 


3 . Fidelity 


Table  II 
hypothesis  tests 
Three  series  of 
hypothesis  tests 
functions  for  Cases 
tests . 


lists  the  results  of  six  series  of 
of  the  types  discussed  in  Section  V-B. 
dispersion  tests  and  three  series  of 
for  the  mean  were  conducted.  Intensity 
III,  IV  and  VI  were  selected  for  these 
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Fcr  each  of  the  three  cases,  a large  number  (500, 
900  or  1000)  of  each  type  of  hypothesis  test  was  performed 
for  a fixed  level  of  significance  a.  The  number  of  times 
each  null  hypothesis  was  rejected  was  recorded  and  the 
statistic  a * {(number  of  re jections)/ ( total  tests)}  was 
computed  for  the  series  of  dispersion  tests  and  the  series 
of  hypothesis  tests  for  the  mean.  If  the  null  hypothesis  is 
true,  then  a is  an  estimate  of  the  significance  level  a. 

From  Table  II  it  appears  that  for  the  hypotheses 
tests  associated  with  Case  III  and  Case  IV,  a gives  a very 
good  estimate  for  a,  the  true  significance  level.  Por  Case 
VI  some  disparities  are  observed.  For  Algorithm  B the  a 
value  for  the  dispersion  test  (i.e.  hypothesis  test  that 
variates  are  Poisson  distributed)  is  0.039  whereas  the 
significance  level  is  a = 0.01.  Algorithm  A yields  a much 
better  a value  of  0.012.  For  the  hypothesis  tests  for  the 
mean,  both  Algorithm  A and  3 yield  a values  almost  twice 
that  of  a • 

These  results,  though  interesting,  are  inconclusive. 
Por  the  dispersion  test,  the  statistic 


is  only  approximately  distributed  as  a chi-sguared  random 

variable.  It  has  been  shown  by  FISHER  [Ref.  4]  that  for 

Poisson  random  variables  with  low  expectations,  hypothesis 

tests  based  on  the  chi-sguared  distribution  may  produce 

spurious  results.  Also,  we  would  expect  that  if  the 

distribution  of  d is  only  approximated  by  the  chi-sguared 
9 P 

distribution,  that  the  relative  error  between  the 

approximate  and  true  distribution  is  greatest  in  the  tails, 
i.e.  at  high  percentile  values.  Case  VI  was  tested  at  an  a 
level  of  0.01.  It  is  also  the  intensity  function  with  the 
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fewest  expected  number  of  events  of  the  three  intensity 
functions  tested.  Partitioning  of  the  interval  into 
segments  produced  some  segments  which  had  a low  expected 
number  of  events  (see  the  right  tail  of  the  density  function 
pictured  in  Figure  12) . The  results  for  Case  71  rather  than 
indicating  any  serious  flaws  in  either  of  the  algorithms 
suggest  further  research  into  finding  better  ways  of 
hypothesis  testing  for  sparse  event  streams  and  highly 
discriminating  levels  of  significance. 

Results  of  the  hypotheses  tests  on  Cases  III  and  17 
indicate  that  both  algorithms  simulate  event  streams  that 
conform  to  the  desired  hon-homogeneous  Poisscn  processes. 


Table  II  - RESULTS  OP  HYPOTHESES  TESTS  FOR  FIDELITY 


C.  GENERAL  OBSERVATIONS 


1 . Programmi nq  Ease 


The  programming  of  Algorithm  A is  simpler  than  the 
programming  of  Algroithm  3.  Because  the  time-scale 
transformation  technique  handles  all  intensity  functions 
alike,  several  different  cases  need  not  be  identified. 
Algorithm  B must  categorize  the  input  intensity  function 
into  one  cf  six  categories  and  then  must  treat  each  category 
in  a unique  way.  A rough  indication  of  this  difference  in 
the  degree  cf  complexity  is  the  number  of  instruction 
statements  required  by  each  program.  The  program 
implementing  Algorithm  A required  142  instructions.  The 
program  for  Algorithm  B required  364  instructions. 

2.  Exact  Method  Versus  Approximate  Method 


Algorithm  B employs  an  exact  method  in  generating 
the  non-hcmogeneous  Poisson  event  stream.  It  is  exact  in 
the  sense  that  all  event  times  are  calculated  directly  and 
the  precision  of  those  calculations  are  limited  only  by  the 
physical  constraints  of  the  computer. 

Algorithm  A employs  a Nevton-Raphson  iterative 
method  tc  approximate  event  times  on  the  interval.  The 
stepping  criterion  for  the  technique  is  limited  by  machine 
precision  considerations.  Also  the  stopping  rule  does  not 
give  a firm  account  of  the  magnitude  of  error  that  can  be 
expected  in  the  event  times.  The  epsilen  criterion  is 
applied  tc  the  value  of  the  integrated  intensity  function 
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and  not  to  the  abscissa,  (or  time  interval  axis) . Por  the 
integrated  intensity  function,  ?(t),  it  is  conjectured  that 
given  a function  value  Pft^)  = u^,  if  a t = t*  can  be  found 
such  that  | F { t • ) - P(t^H  ^ e that  t*  is  very  close  to  t^. 
Although  this  is  a reasonable  assumption  given  the 
well-behaved  nature  of  the  integrated  intensity  function, 
the  problem  of  not  knowing  how  close  t*  is  to  t^  may  tend  to 
reduce  one's  confidence  in  the  algorithm. 


3 • Initialization 

The  initialization  time  required  for  the  programs 
implementing  the  two  algorithms  has  not  yet  been  mentioned. 
Computation  time  comparisons  did  not  include  compilaticn  or 
linkage  times  required  for  the  two  algorithms.  These 
initialization  times  were  significantly  different,  being 
approximately  16  seconds  for  Algorithm  B and  only  8 seconds 
for  Algorithm  A.  Should  simulation  of  only  one  or  two  event 
streams  be  required  it  is  likely  that  Algorithm  A may 
actually  require  less  total  time  than  Algorithm  B.  The 
difference  between  the  initialization  times  could  probably 
be  reduced  if  the  programs  were  to  be  rewritten  in  Assembler 
language . 


D.  CQNCLCSICN  AND  RECOMMENDATIONS 


1 . Conclusion 


For  general  usage,  the  Poisson-decompositio n and  gap 
statistic  method  of  LEWIS  and  SHEDLER  [Ref.  11]  is  clearly 
superior  to  the  time-scale  transformation  method  in 
generating  a non-hoaogeneous  Poisson  process  with  a 
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degree-twc  exponential  polynomial  intensity  function 


The  main  advantages  of  the  Poisson-decomposition  and 
gap  statistic  algorithm  are  its  speed  and  its  qualification 
as  an  exact  method  of  variate  generation. 

The  time-scale  transformation  method  enjoys  some 
advantage  in  core  storage  requirements  and  in  lower  program 
initialization  time.  These  advantages  are  not  sufficient  to 
overcome  the  relative  deficiencies  of  much  greater 
computation  time  and  uncertainty  about  the  accuracy  of  the 
approximating  mechanism  for  determining  event  times. 

2 . Recommendations 

The  full  potential  of  the  Poisson-decomposition  and 
gap  statistic  algorithm  can  not  be  realized  until  it 
includes  a better  method  for  selecting  an  upper  limit  value 
c*  for  Cases  IV,  V and  VI  intensity  functions.  Algorithms 
for  choosing  a c*  value  that  will  be  close  to  the  maximum 
value  for  the  function  x*(t)  should  be  developed  and 
incorporated  into  Algorithm  B.  A single  variable  search 
technique  (such  as  a Golden  Section  search  or  Bisection 
search)  for  finding  the  maximum  value  of  a unimodal  function 
may  be  appropriate. 

The  computer  program  for  Algorithm  B should  be 
rewritten  in  Assembler  language  in  order  to  reduce  program 
initialization  time.  The  program  could  then  be  developed 
into  a library  routine  for  general  use. 

The  question  of  fidelity  of  the  simulated 
ncn-homogenecus  Poisson  process  to  the  true  theoretical 
process  should  be  investigated  further.  Of  special  interest 
is  the  question  of  how  well  the  algorithm  simulates  sparse 
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event  streams.  Methods  for  conducting  the  dispersion  test 
and  the  hypothesis  test  for  the  mean  at  very  high  levels  of 
significance  (i.e.  a - .01,  .005)  for  intervals  with  a low 
mean  number  of  events  are  needed. 
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APPENDIX  A 

PROOF  OF  THE  VALIDITY  OF  SCALING  THE  INVERSE  DISTRIBUTION 

FUNCTION 


Proposition:  Let  T be  a continuous  random  variable  with 

distribution  function  F (•)  , such  that 

T 


FT(t)  = 0 


- _ A (t)  - A (0) 

F*‘  ’ ■ “o 


FT(t)  = 1 


° < t < t0 


t > t. 


Let  Y = ^qF ,p (T)  such  that  T € (0,tQ  ) and  is  a real  number. 

Then  the  random  variable  Y is  uniformly  distributed  cn  the 

interval  (0 , u j . 

0 

Proof : PTC)  maps  the  line  segment  [0,tQ]  onto  the  interval 
[0,1].  If  Ft(«)  is  strictly  increasing  on  (0,tQ)  then 
F”1  (•)  exists  on  the  interval  (0,1)  and  maps  (0,1)  onto 
(0,tQ)  . Nov, 


Fy(y)  = P[Y  < y]  = P[UoFT(T)  < y] 

Z. 

‘0 


= P[F  (T)  < -*->  ] 

1 Uf 


= P [T  < F"1  (-^-)  ] 
“ T y0 

= Ft{F~1 (^-) > 


^0 


Vyl  - £ 


0 < y < u, 
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Therefore, 


dFY(y) 

3y 


fY(y)  = u, 


o < y < u. 


and  Y is  uniformly  distributed  on  the  interval  (0,  Uq)  • 


APPENDIX  B 

EHROH  FUNCTION  AND  DAWSON'S  INTEGRAL  REPHESENTATION  CP  THE 
INTEGRATED  INTENSITY  PUNCTION  A (t) 


The  standard  fora  for  the  error  function  is: 

- t 2 
2 f -u  , 

— / e du  . 

/rr  0 


Dawson's  integral  is  usually  represented  as: 

_ t 2 
t / e u du  . 

0 


Both  of  these  integrals  Bay  be  calculated  to  any  desired 
degree  cf  accuracy  by  using  the  exponential  series 
expansion. 


*X=  Sr 


2 2 
u -u 

Thus,  ® and  e aay  be  written  as 


2 « , 2,n  00  2n 

ou  - ? (u  ) _ 7 y 

‘ n n ! n! 

n=0  n=0 
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T 


and 


e-u2  = l (-u2)n  = l (-l)n(u2n) 


n n! 

n=u  n=u 


respectively.  Integrating  eu  yields 


,t  2 .t  ® 2n 

/ eu  du  = / Y ^ — du 

0 0 n=0  ^ * 

" t 2n 

= l r du 

n=0  0 n! 

oo  ^.2n+l 

ln  (2n+l)  n i 

n=u 

-2 

Multiplying  by  t results  in 


- t 2 °°  ^211-1 

c l e du  = nIo^^ 


and  the  series  representation  for  Dawson's 
revealed.  This  series  will  converge  for  all  t, 
value  of  the  integral  increases  very  rapidly  as 


integral  is 
although  the 
t increases. 
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Using  the  sane  argument  it  is  easily  shown  that  the  series 

. t -u^ 

representation  for  / g e du  is 


" (-1) n t2n+1 
ni0  (2n+l) n! 


The  error  function  is  obtained  by  multiplying  the  integral 
by  the  constant  2//tF  ,i.e.. 


jr  f e-u2 


dU  = -J-  l 

/tt  n=0 


n fc2n+l 

On+UnT 


Although  the  series  expansion  may  be  used  to  calculate  both 
functions,  there  are  very  efficient  computer  library 
routines  available  for  computing  Dawson's  integral  and  the 
error  function. 

The  FORTRAN  supplied  procedures  EP.F  and  DERF  [Ref.  15] 
evaluate  the  above  function  for  input  values  of  t. 

The  IMSL  (International  Mathematical  and  Statistical 
Libraries,  Inc.)  FUNCTION  MMDAW  [Ref.  6]  evaluates  Dawson's 
integral  for  input  values  of  t. 

It  remains  to  be  shown  that  the  intensity  functior 

2 

X (t)  = exp{ot^  + a^t  + a^t  } may  be  integrated  over  the 
interval  (0,tQ]  using  either  ERF  (or  DERF)  or  MMDA'4. 


I 
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M. 


Consider , 

fc0  fc0  - 
A(tn)-A(0)  * / X(t)dt  * / exp(an  + a. t + a_t  ;dt 

o o u 4 


By  completing  the  square  of  the  exponent  the  expression 
beccaes 

A<V  * l ° exp  { “O  ' 4^}  eXP  { °2(t  + 2^’}  dt 


( A (°)  = since  0 is  the  left  end  point  of  the  interval 
over  which  \ (t)  is  defined)  . 


For  a, > 0 , 


A ( tQ ) = exp  < aQ 


Letting, 


/ 0 exp  | -SJ  (t  ♦ dt 


u ‘ ^2  (t  + SJ1 


du  = /a2  dt 


K'  = exp  / a.  - 


0 4a- 


and  adjusting  the  limits  of  integration,  the  expression 
becomes: 


A(t0)  = 


b 2 


/ eu  du  = 


'a2  a 


b 2 a 2 ) 

/ eu  du  - / eu  du  [ 
0 o ) 
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where. 


and 


a = 


b = 


W 


The  expression  atove  can  be  rewritten  as  follows: 


A(t0)  = 


_K_ 

'/cT 


{b2b~2  / eu2du  - a2a~2  / eu  du> 
0 0 


Inside  the  brackets  are  two  Dawson  integral  expressions 
multiplied  by  constants.  Outside  the  brackets  is  just  one 
more  easily  determined  constant.  Using  MIJDA»  twice  and 
multiplying  by  the  appropriate  constants  will  give  the 
desired  integral  value. 


For  <*2  < 0, 

A(tQ)  = exp  | aQ 


/ exp 
0 


/ exp 
0 


al  2 


" [/| a2 ' 


(t  + 


Substituting  as  before  gives 


A ( tQ ) 


K 

/R“T 


du 


du  - 
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and  the  error  function  can  be  recognized  within  the 
brackets. 

Two  error  function  calculations,  a subtraction  and  a 

multiplication  by  a constant  are  sufficient  to  determine 

A(t  J . 
o 
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APPENDIX  C 


COMPARISON  OP  GAP  STATISTICS  ALGORITHM  AND  TIME-SCALE 
TRANSFORMATION  ALGORITHM  POR  SIMULATION  OF  NON- HOMOGENEOUS 
POISSON  PROCESSES  WITH  LOG-LINEAR  INTENSITY  FUNCTIONS 


LEWIS  and  SHEDLER  [Ref.  13]  present  two  algorithas 
for  simulation  of  non-homogeneous  Poisson  processes  with 
intensity  function  X (t)  = exp(BQ  *6^)  . The  first  algorithm, 
which  uses  a time-scale  transformation  is  easy  to  employ 
since  the  inverse  of  the  integrated  intensity  function  is 
known  explicity.  The  second  algorithm  uses  gap  statistics 
to  generate  event  streams. 

Ecth  methods  are  exact.  That  is,  other  than  the 
physical  precision  constraints  inherent  to  the  computer,  no 
approximations  are  necessary  in  generating  event  times. 
However  the  gap  statistics  algorithm  obviates  the  need  for 
ordering  an  array  of  variates.  Also,  the  gap  statistics 
algorithm  takes  advantage  of  a fast  standard  exponential 
variate  generator,  (Random  Number  Generator  Package 
LLRANDOM,  Ref.  12),  thereby  avoiding  the  time  consuming 
computation  cf  taking  logarithms.  The  time-scale  algorithm 
must  perform  both  the  ordering  of  variates  and  the 
computation  cf  logarithms. 

Beth  algorithms  were  programmed  and  tested  for 
computational  speed  efficiency.  The  gap  statistics 
algorithm  was  roughly  twice  as  fast  as  the  time-scale 
transf or maticn  algorithm. 
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The  SOBEOOTINE  NHPP2  uses  the  gap  statistics 
algorithm  within  the  Poisson-decomposition  and  gap  statistic 
algorithm  (Algorithm  B)  . It  is  the  presence  of  SOBRCOTINE 
NHPP2  which  markedly  reduces  the  ordering  requirements  of 
Algorithm  B from  those  necessary  in  Algorithm  A. 
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APPENDIX  D 


POISSON  VARIATE  GENERATION 


A.  BACKGROUND 

Acknowledgement:  The  content  of  this  appendix  is  a 
paraphrase  cf  portions  of  Chapter  11  of  an  unpublished  book 
by  J.  H.  AHRENS  and  0.  DIETER  [Ref.  1]. 


The  Pcisson  distribution  has  the  probability  mass 
function  , . 


where  is  the  probability  that  exactly  i events  are 
observed  to  cccur  in  a unit  time  interval,  given  that  events 
arrive  at  an  average  rate  \ . 

The  intervals  between  events  in  a Poisson  process  of 
rate  \ are  independent  and  have  an  exponential  distribution 
with  mean  \/\,  i.e. 

f(x)  = Ae 

The  probability  integral  transform  method  can  be  used  to 
obtain  a sample  of  exponential  variates  from  a sample  of 
uniform  (0,1)  variates,  u ,u  ,.  • • • since 

X1  = ~ X Ul'  x2  = ” x u2  ' ' * * 


A X rate  event  stream  can  then  be  simulated  from  the 
exponential  variates  using  the  following  interpretation: 


1st  event  occurs  at  x 


2nd  event  occurs  at  x ♦ x 

1 2 

etc . 


This  property  yields  two  simple  methods  for  generating  a 
realization  cn  the  Poisson  random  variable.  These  are  the 
additive  method  and  the  multiplicative  method. 


B.  THE  ADDITIVE  HETHOD  ?0R  GENERATING  A POISSON  VARIATE 


The  number  k of  units  arriving  in  a unit  interval  in  a 
Poisson  process  with  rate  X must  be  found.  This  will  be 
the  k for  which 


X1  + X2  + * * ’+  £.  1 


*1  + x2  +'”+  xk  + xk+l  =•  1 


or  equivalently,  for  which 


-Jin  - Jin  U2  -•••-  Jin  £ X 


-Jin  - Jin  u2  -•••-  Jin  - £n  uk+1  >_  * (1) 


(Note:  Since  u^  < 1 Vi  , -‘■In  u^  > 0) 


Thus  by  using  uniform  variates,  computing  logarithms, 
adding  and  counting,  the  Poisson  variate  is  obtained.  The 
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problem  with  the  additive  method  is  that  if  x is  very 
large,  say  2000,  it  requires  considerable  effort  to  generate 
a single  Poisson  variate.  Although  addition  is  extremely 
fast  on  a computer,  the  computing  of  logarithms  is 
relatively  slow. 


1 ' 

* 1 
3 1 


C.  THE  MULTIPLICATIVE  METHOD  FOR  GENERATING  A POISSON 
VARIATE 


Multiplying  (1)  by  -1  gives: 

in(u1*u2 uk) 

and 

*n(uru2 uk*uk+i) 

or  equivalently 


> -X 


< -X 


and 


u^*u2 e 


-X 


VU2 


V'W  6 


So  by  generating  uniforms,  successive  multiplication  and 
counting,  the  Poisson  variate  k is  produced.  The 

multiplicative  method  eliminates  the  need  for  computing 
logarithms.  However  trouble  is  encountered  for  large 
values  since  e”^  becomes  very  small  as  X increases. 
Underflow  problems  occur  for  X values  of  approximately  175 
and  larger  on  the  IBM  360/67  computer  system  using  single 
precision  word  lengths.  The  precision  problem  may  be 

overcome  by  partitioning  the  Poisson  random  variable  into 
the  sum  of  two  or  more  Poisson  random  variables.  However 
the  average  number  of  multiplications  and  uniform  variates 
required  is  still  proportional  to  X . More  efficient 


I 


96 


methods  exist  such  as  the  Gamma  method 


0.  THE  GAMMA  METHOD  FOR  GENERATING  A POISSON  VARIATE 

(The  following  discussion  is  a condensed  account  of  Ahrens 
and  Dieter's  presentation;  Ref.  1,  pp. 11-20,21,22) 

In  order  to  obtain  a sample  k from  the  Pcisson 
distribution  of  mean  y , select  a positive  integer  n 
(typically  n is  a little  smaller  than  M).  Then  take  a 
sample  x from  the  Gamma  distribution  of  parameter  n. 

(1)  If  x > y return  a 

sample  k from  the  binomial 
distribution  with  parameters 
n-1  , V /x 

(2)  If  x S U take  a sample 

j from  the  Poisson 

distribution  of  mean  y-x  and 
return  k «-  n ♦ j. 


The  sample  x simulates  the  n-th  event  (arrival)  in  a 
Poisson  process  of  rate  1.  If  x > n then  there  are  n-1 
arrivals  in  the  interval  (0rx) , and  each  of  these  has  a 
probability  cf  y/x  of  being  below  m [Case  (1)  ].  If 
x S y then  the  n simulated  arrivals  are  all  before  U, 
and  the  sample  j indicates  the  additional  events  between  x 
and  y [Case  (2)  ]. 

(At  this  point  Ref.  1 includes  a formal  proof  of  the 
procedure,  which  will  be  omitted.) 


The  explicit  algorithm  relies  on  an  efficient  method  for 
saapling  free  the  Gaaea  distribution.  (Such  a aethod  has 
been  iapleaented  as  9.  H.  CHURCH  Coaputer  Center  library 
routine  GAHA  by  D.  9.  Robinson  and  P.  A. 9.  Levis  and  is 
docuaented  in  Reference  10.)  The  constant  n is 
arbitrarily  chosen  as  n * [0.875  y ]•  This  avoids  the  Case 
(1)  alaost  ccapletely  if  y is  large  and  a simple  aethod  for 
saapling  froa  the  binoaial  distribution  becomes 
sufficient:  the  Bernoulli  process  is  siaulated  in  Steps 
8-10  belcv.  The  algorithm  for  sampling  froa  the  Poisson 
distribution  with  mean  y - x in  Case  (2)  is  the  simple 
multiplicative  method  explained  above  and  is  employed  in 
Steps  3-5.  However,  if  y - x is  still  large  (larger  than 
the  "cut-cff  point”  c)  the  entire  procedure  is  iterated  with 
a new  auxiliary  sample  x from  a Gamma  distribution  of  mean 
n'  * [0.875  (y  - x) ] etc. 


Algorithm  - Jhe  Gamma  Method,  Ahrens  and  Dieter 


1.  Initialize  k «-  0,  w «-  y. 

2.  If  w £ c go  to  6 (c  = 16  was  determined 
empirically  by  Ahrens  and  Dieter  to  be 
reasonable) . 

3.  (Start  Case  (2)).  Set  p •<-  1 and  calculate 

. -w 
b *■  € . 

4.  Generate  u and  set  p «-  p»u.  If  p<b  deliver 
k. 
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Increase  k «-  k+1  and  go  to  4 


6.  Set  n «-[dw]  where  d = 7/8  (d  value  was  also 
selected  empirically) . Take  a sample  z from  the 
Gamma  (n,1)  distribution.  If  x > w go  to  8. 


7.  Set  k = k + n,  w «-  w-x  and  go  to  2. 

8.  (Start  Case  (1)).  Set  p +-  w/x. 

9.  Generate  u.  If  u<  p increase  k k ♦ 1. 

10.  Set  n «-  n-1.  If  n > 1 go  to  9. 

11.  Deliver  k. 

Ahrens  and  Dieter  claim  that  the  computation  time  for 
the  Gamma  method  grows  ultimately  like  a + 8 In  y . 
Computation  time  for  the  additive  and  multiplicative  methods 
grew  like  u . Therefore  for  the  large  u values  typically 
demanded  in  the  simulation  of  non-homogeneous  Pcisson 
processes,  the  Gamma  method  is  signally  superior  in  terms  of 
computation  speed. 


APPENDIX  E 


GBAPHICAL  PRESENTATION  OP  NEWTON-R APHSON  METHOD 


/ 

A.  GENERAL  PRINCIPLES 


Figure  E-l 

Consider  the  function  F (t)  pictured  above  in  Figure  E-1. 

The  objective  is  to  find,  for  a known  value  u^  on  the 

vertical  axis,  a unique  corresponding  value  t^  such  that 

P (t^)  * u^.  If  an  explicit  expression  for  the  inverse  of 

P(«)  exists,  the  calculation  nay  be  nade  directly,  i.e., 

t * F“l(u  ).  Otherwise,  nuaerical  aethods  aust  be  used  to 
i . i 
appxoxiaate  t . 

i 
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; 


6i(t)  * F(t)  - ui 


Form  the  new  function  G^{t  ) - F(t  ) - (see  Figure 

E-2)  . In  effect,  the  horizontal  axis  has  been  displaced 
upward  a distance  u^  Now  if  a t*  can  be  found  such  that 
G^(t*)  = 0,  then  F(t*)  = u^  and  t*  = tif  the  desired  value. 


Assume  that  an  initial  approximation  for  t*  can  be  made, 
say  tj  . If  (t?  ) and  G^  (t!  \ = gi  (t^  ) can  be  computed 
we  can  write 


g . (t ! ) 


{TT 


(tf) 


where  is  the  intercept  of  the  line  tangent  to  G^  (•)  at 

Gi  ) with  the  t-axis.  Solving  the  above  expression  for 


t* 

2 


gives: 


*2 


G (tp 

t ' - ~ 1 

1 


It  is  evident  from  the  graph  that  t^ 
to  the  unknown  t*. 


is  closer  than  t* 
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Repeating  the  procedure  (Figure  2-3)  , using  t£  as  a new 
approximation  will  yield  t^  a still  better  approximation  of 


t*. 


In  general, 
a clo 

the  relationship 


given  an  approximation  t£ 


obtained,  a closer  approximation  t£+1 


has  been 
can  be  calculated  by 


’ti+I  ~ tk“ 


Gi(V 
1 V 


It  can  be  shown  that  successive  approximations  will  converge 
to  the  value  t*  provided  G^(t)  satisfies  certain  criteria 
[Ref.  5,  F.  449, 450  ]. 

Since  G^  (t)  has  the  form  of  a distribution  function  F (t) 
it  satisfies  all  the  conditions  for  convergence  except  for 
the  case  where  the  interval  [t^  , t/J  contains  an  inflection 
point.  If  the  inflection  point(s)  can  be  identified,  this 
troublesome  situation  can  easily  be  avoided. by  handling  the 
function  G^  (•)  in  discrete  sections  over  the  discrete 
intervals,  <0,t«  ),  (t*«  ,t«  ) , . . . , ) where  the  t •» 
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mm 
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(j  * 1,  . . . , k)  are  values  such  that  Gi(t^  ) are  all  the 

inflection  points  of  G^(«)  over  (0,tg].  (Por  the  special 
case  under  consideration  at  nost  one  inflection  point  will 
be  encountered.  It  is  the  maximum  or  ainiaun  of  the 
intensity  function.) 

Successive  approxinat ions  are  computed  until 
| G i (t  ) | < e at  which  tine  t£  is  assumed  to  be  close 
enough  to  t*. 

B.  INITIAL  APPFOXIN ATIONS 

The  problem  of  obtaining  a good  initial  approximation 
for  each  t^  was  solved  in  the  following  manner  (see  Figure 
E-4  below)  : I 


Figure  E-4 
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The  interval  (0,t(j]  was  divided  into  several  segments 
(0#t  ],  (x;L»t2]/  etc.  Then  the  values  F(Xj)  (j  = 1,2,...) 

were  computed.  Bach  resulting  abscissa  and  ordinate  value 
was  stored  in  an  array,  i. e. , 


t Fit) 


0 

0 

T1 

Fftj) 

T2 

F ( Tg ) 

• 

• 

t0 

F(t0) 

For  each  value  u the  array  was  searched  until  two  function 
values  were  located  such  that  F(tj)  < u < F<Tj+i)  • Either 
x j or  Tj+1  was  then  selected  as  the  initial  approximation 
for  the  value  t.  (i.e.  t-  ,)  such  that  F(tJ  * u • . The 

X X f X XX 

function  G^(t  ) = F(t  ) - u^  was  then  fcried  and  the 
Newton-Raphscn  iterations  performed  until  the  epsilon 
criterion  was  satisfied. 
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SUBROUTINE  MTOINV 
PURPOSE 

SIMULATES  A NON-HOMOGEN EOUS  POISSON  PROCESS  WITH 
QUADRATIC  EXPONENTIAL  INTENSITY  FUNCTION  OVER  A 
GIVEN  INTERVAL  USING  A TIME-SCALE  TRANSFORMATION 
OF  A HOMOGENEOUS  POISSON  PROCESS. 


CALL  MTOINV  < A, A1 , A2, EL t ER, I S t I I t M, I ER ) 

DESCRIPTION  OF  PARAMETERS 

A - CONSTANT  IN  INTENSITY  FUNCTION. 

Al  - 1ST  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

A2  - 2ND  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

EL  - LEFT  END  POINT  OF  INTERVAL. 

ER  - RIGHT  END  POINT  OF  INTERVAL. 

IS  - RANDOM  NUMBER  SEED.  ANY  INTEGER  WITH  NINE 
OR  LESS  DIGITS. 

II  - 0 FOR  TIMES  OF  EVENTS. 

I FOR  TIMES  BETWEEN  EVENTS. 

M - RETURNED  VALUE.  EQUALS  NUMBER  OF  EVENTS  IN 
NON-HOMOGEN EOUS  POISSON  PROCESS. 

I ER  - ERROR  CODE.  MAY  BE  WRITTEN  ON  DEMAND.  I ER 
EQUALS  THE  NUMBER  OF  TIMES  PROGRAM  WAS 
FORCED  TO  ABORT  SIMULATION  AND  START  AGAIN 
DUE  TO  EXCESSIVE  EVENTS  (MORE  THAN  5000). 


CCMMENT  S 

CALLING  PROGRAM  MUST  HAVE  A COMMON  REGION,  BLOKl, 
OF  DIMENSION  (5000). 

EXAMPLE:  DIMENSION  T ( 5000 ) 

COMMON/ BL0K1/T 

CALLING  PROGRAM  MUST  USE  THE  FOLLOWING  JCL  CARDS 

//  EXEC  FORTCLG,  IMSL*DP 
//FORT. SYS  IN  DO  * 

TIMES  TO  EVENTS  OR  TIMES  BETWEEN  EVENTS  WILL  BE 
STORED  IN  CELLS  Til)  THROUGH  TIM). 


SUBROUTINE  MTDINV  ( A, Al ,A2,EL, ER, IS, 1 1,  M, IER) 
DIMENSION  TAU15000) , 8RKT<5000,2) 

COMMON  /8L0K1 / TAU/BLQK2/BRKT 

CALL  OV  FLOW 

IER  = - 1 

IER  * IER+1 

BRKT( 1,1)  = 0. 

BRKTI  1,  2)  * 0. 

INTEGRATE  INTENSITY  FUNCTION 

CALL  HELP  ( A, Ai,A2,EL,ER,PAR) 

PARI  * PAR 

GENERATE  POISSON  VARIATE  M 

CALL  PVAR  ( IS,  PAR  1 * M ) 

IF  ( M. GT. 5000)  GO  TO  1 

SET  EPSILON  VALUE  FOR  NEWTON-RAPHSON  STOPPING  RULE 
EPS  * FLOAT(M+l)*C. 00001 
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GENERATE  M UNIFORM  VARIATES;  OROER  AND  SCALE 

CALL  RANOOM  (IS,TAU,M) 

CALL  PX  SORT  <TAU,l,M) 

DO  2 1*1, M 

TAU< I ) * T AU( I ) *PAR 

2 CONTINUE 

DETERMINE  INTERVAL  PART  ION  I NG  SCHEME 
FOR  SUBROUTINE  BNCHMK 

12  * M/4 

OIV  * AMAX0I10, 12  ) 

IDIV  = INT( DI V+2 ) 

DELTA  = (ER-ED/DIV 
ENDL  * EL 

COMPUTE  ARRAY  OF  ABSCISSAS  AND  ORDINATES  FOR 
INTEGRATED  INTENSITY  FUNCTION 

CALL  BNCHMK  ( A, AI , A2 , ENDL, DELTA , IDIV ) 

BEGIN  NEWTON-RAPHSON  ITERATIONS 

K * 0 
J = 1 

3 K = K+l 

4 IF  (TAU(K).LT.BRKT( J ,2) ) GO  TO  5 
J = J+l 

GO  TO  4 

5 IF  (A8S<BRKT<J,2)-TAU(K>).LT.ABSl9RKT(J-i,2)-TAU(K>)) 
*GG  TO  6 

T = BRK  T(  J-  l » 1 ) 

A CD  * ( BRKT(J-1,2)-TAU(K) ) / EVAL ( A ,A 1 ,A 2 , T ) 

GO  TO  7 

6 T * BRKT(  J , 1) 

ADD  * ( BRKT(J,2)-TAU(K) )/EVAL(A, A1,A2,T) 

7 IF  ( ABS ( ADD ) . LT. E PS)  GO  TO  8 
T = T-ADO 

CALL  HELP  ( A, AI,A2, EL,T,X) 

ADO  * ( X-TAU(K)  )/EVAL(A,Al  ,A2,T) 

GO  TO  7 

8 TAU(K)  = T 

IF  (K.EQ.M)  GO  TO  9 
GO  TO  3 

9 IF  (II.EQ.O)  RETURN 

TIMES  BETWEEN  EVENTS  ARE  REQUESTEO-CALCULATE  THEM 

S = TAU(l) 

TAU(l)  * TAU<1)-EL 

DO  10  1=2, M 
SI  * TAU( I ) 

T AU(  I ) * TAU(  I)-S 
S * SI 
10  CONTINUE 

RETURN 

ENO 


SUBROUTINE  HELP  CALCULATES  THE  VALUE  OF  THE 
INTEGRATED  INTENSITY  FUNCTION  OVER  ANY  INTERVAL  (0,T) 
OVER  WHICH  THE  NHPP  OCCURS 

SUBROUTINE  HELP  ( A, Al , A2, EL, ER, XX ) 

DOUBLE  PRECISION  MMD AW , BB » A A 
Z = SORT ( ABS( A2 ) ) 

Y * ( Al*Z) /( 2.*A2 ) 
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* A-A1*A1/(4.*A2) 

■ EXP(CC  )/Z 
(A2.LT.0.)  60  TO  2 

* DE  XP( A A**2 ) * HMD AW ( AA) 


OC  1 
BINC 
8B  * Z*BINC+Y 

Q2  = DEXP(BB**2)*MMDAW(B8) 
XX  = CC*(  Q2-Q1 ) 

BRKT ( 1,1 ) * BINC 
BRKT(  1,2)  = XX 
1 CONTINUE 


1 = 2 » 1 0 IV 
* BINC+OELTA 


GO  TO  3 

2 CC  = CC*. 8862269 


00  3 I *2, 10 IV 
BINC  * BINC+OELTA 
BB  = Z*B INC +Y 

XX  * CC*( OERF ( BB )-OERF( AA) ) 
BRKT (1,1)  = BINC 
BRKT ( 1,2)  * XX 
3 CONTINUE 


RETURN 

ENO 


FUNCTION  EVAL  COMPUTES  THE  VALUE  OF  THE 
INTENSITY  FUNCTION 

FUNCTION  EVAL  (A,Al,A2,T) 

EVAL  = A+A1*T+A2*T**2 
EVAL  * EXP  (EVAL) 

RETURN 

ENO 
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SUBROUTINE  OEGTWO 


r 


PURPOSE 

SIMULATES  A NCN-HOMOGENEOUS  POISSON  PROCESS  WITH 
QUADRATIC  EXPONENTIAL  INTENSITY  FUNCTION  OVER  A 
GIVEN  INTERVAL  USING  THE  POISSON-DECOMPOSITION 
AND  GAP  STATISTIC  ALGORITHM. 

USAGE 

CALL 


DEGTWO  ( I S,A,A1 ,A2,EL,ER,II ,Nt I ER  > 


DESCRIPTION  OF  PARAMETERS 


ANY  INTEGER  WITH  NINE 


IS  - RANDOM  NUMBER  SEED. 

OR  LESS  DIGITS. 

A - CONSTANT  IN  INTENSITY  FUNCTION. 

A1  - 1ST  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

A2  - 2ND  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

EL  - LEFT  END  POINT  OF  INTERVAL. 

ER  - RIGHT  END  POINT  OF  INTERVAL. 

II  - 0 FOR  TIMES  OF  EVENTS. 

1 FOR  TIMES  BETWEEN  EVENTS. 

N - A VECTOR  OF  LENGTH  5.  N(l> 

CONTAIN  NUMBERS  OF  EVENTS 
COMPONENTS  OF  THE  DECOMPOSED  INTENSITY 
FUNCTION.  N ( 5 ) CONTAINS  THE  TOTAL  NUMBER 
OF  EVENTS  IN  THE  NON-HOMOGENE OUS  POISSON 
PROCESS . 


THROUGH  N<4) 
FROM  VARIOUS 


COMMENTS 

CALLING  PROGRAM  MUST  HAVE  A COMMON  REGION»  HOLD, 
OF  DIMENSION  (5000),  AND  AN  INTEGER  ARRAY  OF 
DIMENSION  (5). 

EXAMPLE:  DIMENSION  TI5000),  N(5> 

COMMON/HOLD/T 

CALLING  PROGRAM  MUST  CONTAIN  THE  FOLLOWING 
ASSIGNMENT  STATEMENT: 

M=N(5> 

CALLING  PROGRAM  MUST  USE  THE  FOLLOWING  JCL  CARDS 

//  EXEC  FORTCLG, I MSL=DP 
//FORT.SYSIN  DD  * 

TIMES  TO  EVENTS  OR  TIMES  BETWEEN  EVENTS  WILL  BE 
STORED  IN  CELLS  T(l)  THROUGH  TIM). 


SUBROUTINE  DESTWO  (I S , A , A1 , A2 , EL , ER, I I , N , I ER) 
DIMENSION  T IMES  ( 5 000  ) , T(5000),  N<5),  P<5) 
COMMON  /MIKE/  TIMES/HOLD/T 
CALL  OVFLOW 

INITIALIZE  VARIABLES 


P(l) 
P ( 2 ) 
P ( 3 ) 
P<4) 
P ( 5 ) 


A 

A1 

A2 

0. 

0. 


DO  1 1*1,5 
1 N( I)  = 0 

IF  RATE  FUNCTION  IS  LESS  THAN  DEGREE  TWO, 
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USE  NHPP2  ROUTINE  ONLY 

IF  (A2.EQ.0.)  GO  TO  2 
GO  TO  4 

2 CALL  NHPP2  ( I S, EL » ER, A, A1 » 1 1 ,N1 , 1 ER ) 

N<  5 ) = N 1 

IF  (N1.EQ.0)  RETURN 

DO  3 I*1»N1 
TIMES ( II  * T{ I) 

3 CONTINUE 
RETURN 

DETERMINE  COEFFICIENTS  FOR  MODIFIED 
DEGREE  ONE  RATE  FUNCTION 

4 TEST  = -Al/(2.*A2) 

TINT  = ER-EL 

IF  (AI. GE. 0.. ANO. A2. GT.O.)  GO  TO  5 
GC  TO  6 

58=  A-A2*T I NT**2 
81  * A1+2.*A2*TINT 
GO  TO  10 

6 8 * A 

IF  ( ( A I .LE.O. • AND .A  2. LT. 0. ) . OR . ( AI .GT.O . . AND. A2 .LT .0 . 
*. ANO. TEST. GE. TINT  ) ) GO  TO  7 
GO  TO  8 

7 81=*  A L +A  2*TINT 
GO  TO  10 

8 IF  (AI. GT. 0. . ANO. A2.LT.0.. AND.TEST.LT. TINT)  GO  TO  9 
81  - AL 

GO  TO  10 

9 81  = Al/2. 

MUST  THE  INTERVAL  0E  PARTITIONED? 

10  IF  (Al*A2.LT.O.. ANO.TEST.LT. TINT)  GO  TO  11 
ERNEW  = ER 

GO  TO  12 

11  ERNEW  * TEST+EL 

GENERATE  OEGREE  ONE  NHPP  ON  INTERVAL 

12  88  * 8 

?B1  s g ] 

ALL  NHPP 2 ( I S » E L * ERNEW, BB, BBl ,0»N1, IER) 

NIL)  * Nl 

IF  ( N ( 1 ) . EQ.O)  GO  TO  14 


T I MES ( I)  = T(  I) 
13  CONTINUE 


COMPUTE  LENGTH  OF  INTERVAL  AND  DETERMINE  VALUE 
OF  CSTAR  FOR  USE  IN  REJECTION  ROUTINE 

14  Q * ERNEW-EL 
El  = A 

E2  * A2*Q**2 
E3  » AI *Q 

E4  = Al**2/(4.*A2  ) 

E5  = Al**2/( 2.*A2  ) 

IF  (AI  .GE.O.. ANO. A2. GT.O. ) GO  TO  15 

IF  (AI. LT. 0. . AND. A2. GT.O. .AND. TEST. GE. TINT)  GO  TO  16 
IF  (Al.LT.O..  AND.  A2.  GT.O.  .ANO.TEST.LT.  TINT)  GO  TO  17 
IF  (AI. LE.O.. ANO. A2. LT. 0. ) GO  TO  18 

IF  (Al.GT.O.. AND. A2.LT.0. .AND. TEST. GE.T INT)  GO  TO  15 
CSTAR  = EXP(E1-E4)-EXP(E1) 

GO  TO  20 

15  CSTAR  * EXP(E1)-EXP(E1-E2) 

GO  TO  20 


110 


nnnnnnn  nono  o ooo  non  ooo  non  o non  ooo  ooo  noon 


16  CSTAR  * EXP(E1+E2+E3)-EXP(E1+E3) 

GO  TO  20 

17  CSTAR  * EXP( El -E4 ) -EX  P( EI-E5) 

GO  TO  20 

18  CSTAR  = EXP ( El )-EXP( E1  + E3+E  2) 

GO  TO  20 

19  CSTAR  * EXP(E1+E3+E2)-EXP(EI) 

COMPUTE  INTEGRAL  OF  MODIFIED  DEGREE  TWO  RATE  FUNCTION 
OVER  INTERVAL 

20  CALL  HELP  (A«A1*A2»EL,ERNEW.PMTR) 

PMTR  = PMTR-(EXP(Bj*(EXP(Bl*ERNEW)-EXPtBl*EL)  ) ) /B  1 

I0ENTIFY  AS  FIRST  SUBINTERVAL 
NOTE  * 1 

GENERATE  REALIZATION  ON  POISSON  (PMTR)  VARIATE 

21  CALL  PVAR  (IS, PMTR, M) 

IF  ( NOTE.EQ.l)  GO  TO  22 
GO  TO  25 

REJECTION  ROUTINE  USED  ON  FIRST  SUBINTERVAL 

22  N ( 2)  = M 
P<4)  * B 
P (5  ) = B1 

IF  (N(2).EQ.O)  GO  TO  24 

CALL  REJECT  ( IS , EL ,C STAR, P , Q, N ( 2) ) 

00  23  I-l.M 
T IMES  ( N ( 1 1 + 1 ) = T ( I) 

23  CONTINUE 

HAS  THE  INTERVAL  BEEN  PARTITIONED? 

24  IF  (ERNEW.EQ.ER)  GO  TO  34 
GO  TO  27 

USE  REJECTION  ROUTINE  ON  SECOND  PART  OF  INTERVAL 

25  N ( 4)  * M 
P (4  ) = B 
P ( 51  a Bi 

IF  NO  EVENTS  OCCURRED  BYPASS  REJECTION  ROUTINE 

IF  (N(4).EQ.O)  GO  TO  35 
Q = ER-ELNEW 

CALL  REJECT  (I  S , ELNE W ,CST AR ,P , Q , N( 4)) 

COPY  TIMES  OF  EVENTS  INTO  'TIMES'  ARRAY 

N4  » NC l»+N(2)+N( 3) 

00  26  I * 1 • M 
T IMES (N4+ I ) * T< I) 

26  CONTINUE 

GENERATION  OF  VARIATES  COMPLETE. 

GO  TO  ORDERING  ROUTINE 

GO  TO  35 

INTERVAL  PARTITION  WAS  REQUIRED.  MUST  NOW 
CONSIDER  SECOND  SUBINTERVAL 

DETERMINE  COEFFICIENTS  FOR  MODIFIED 
DEGREE  ONE  RATE  FUNCTION 
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27  IF  (A1.GT.0..AND.A2.LT.0.)  GO  TO  28 
B * A-A2*TINT**2 

B 1 » A1+2.*A2*TINT 
GO  TO  29 

28  B - A+( Ai/2.)*TINT 
B1  * Ai/2.  + A2*TINT 

29  ELNEW  * ERNEW 

GENERATE  OEGREE  ONE  NHPP  ON  INTERVAL 

BB  » B 
BB1  * Bi 

CALL  NHPP 2 ( IS,ELNEW,ERtBBt8BIt0tN3,IER) 

N (3  ) » N3 

IF  (NO).EQ.O)  GO  TO  31 
N3  » N( 1)+N(2) 

TRANSFER  TIMES  BETWEEN  ARRAYS 

00  30  I»1.N3 
T IMES( N 3+ 1 ) * T( I ) 

30  CONTINUE 

31  Q » TINT 

DETERMINE  VALUE  OF  CSTAR  FOR  USE  IN 
THE  REJECTION  ROUTINE 

E 2 = A2*Q**2 
E3  » A1*Q 

IF  (A1.GT.0..AND.A2.LT.0.)  GO  TO  32 
CSTAR  = EXP(Ei-E4)-EXP(El-E5-E3-E2) 

GO  TO  33 

32  CSTAR  - EX  P( E1-E4 }-EX  P( E1+E3+E2) 

COMPUTE  INTEGRAL  OF  MODIFIED  DEGREE  TWO  RATE 
FUNCTION  OVER  SECOND  INTERVAAL 

33  CALL  HELP  ( A,  Al , A2 , ELNEW,  ER  tPMTR  ) 

PMTR  * PMTR-(EXP<B)*tEXP(Bl*ER)-EXP<Bl*ELNEW)))/Bl 

IDENTIFY  AS  SECOND  SUBINTERVAL 

NOTE  » 2 
GO  TO  21 

PARTITION  OF  INTERVAL  NOT  REQUIRED.  COMPUTE  TOTAL 
EVENTS  AND  SUPERPOSE  TWO  EVENT  STREAMS 

34  N(  5 ) « N<l)+N<2) 

IF  INI  2 ) . EQ.O)  GO  TO  38 
LBGN  » N ( 1 ) +1 
JBGN  * 1 

CALL  COLATE  <LBGN,N<5),1) 

GO  TO  38 

PARTITION  WAS  REQUIRED.  DETERMINE 
AMOUNT  OF  SORTING  NEEOED 

35  N ( 5 ) » N(l)+N(2)+N(3)+N(4)* 

IF  (N(2). EQ.O. AND. N<4). EQ.O)  GO  TO  38 
IF  (NIA).EQ.O)  GO  TO  36 
IF  (N(Z).EQ.O)  GO  TO  37 

MUST  SUPERPOSE  FOUR  EVENT  STREAMS 

LBGN  » NdH-l 
LF  IN  « N(  1 )+N( 2 ) 

CALL  COLATE  ( LBGN ,LFIN, 1) 

LBGN  = LF  I N+N(  3 ) •*•  l 
JBGN  * LFIN+l 

CALL  COLATE  <LBGN,N<5 ) , JBGN ) 
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GO  TO  38 

MLST  SUPERPOSE  FIRST  HALF  OF  ARRAY  ONLY 

36  N2  * N(l)+N(2) 

LBGN  * N( 1 ) +1 

CALL  COLATE  (LBGN,N2,l> 

GO  TO  38 

MUST  SUPERPOSE  SECOND  HALF  OF  ARRAY  ONLY 

37  KK  * N(  l)*Nt2m 

LBGN  = N(l)+N(2)+N(3)+l 
LFIN  = N{ 5) 

CALL  COLATE  (LBGN ,N( 5 > , KK ) 

GO  TO  3 8 

ARE  TIMES  OF  EVENTS  OR  TIMES  BETWEEN  EVENTS  REQUESTED? 

38  IF  (II.EQ.O)  RETURN 

CALCULATE  TIMES  BETWEEN  EVENTS 

S » TIMES  ( 1 ) 

TIMES!  I ) * TI  ME  S(  U-EL 
N5  » N( 5 ) 

00  39  I-2.N5 
SI  * TIMES(I) 

TIMES(I)  » TIMES( I)-S 
S = SI 

39  CONTINUE 
RETURN 
END 


SUBROUTINE  NHPP2  SIMULATES  A NON-HOMOGENEOUS 
PCISSON  PROCESS  WITH  A LOG-LINEAR  INTENSITY 
(RATE)  FUNCTION 

SUBROUTINE  NHPP2  ( I S , EL r ER r A, Al f 1 1 , N, I ER ) 

DIMENSION  T (5000 ) 

COMMON  /HOLO/  T 

CALL  OVFLOW 

INITIALIZE  VARIABLES 

IER  « 0 
TINT  a ER-EL 
A » EXP ( A+AI*EL ) 

IS  THE  POISSON  PROCESS  HOMOGENEOUS? 

IF  (AL.EQ.O.)  GO  TO  3 

PAR  * ( A*< EXP(TINT*AI )-l. ) )/Al 

IF  (AI.GT.O.)  30  TO  1 

I FLAG  = 3 

GO  TO  2 

A * a*EXP(  TINT*AI ) 

Al  = -Al 
I FLAG  a 2 

COMPUTE  PARAMETERS  OF  BOTH  POISSON  RANDOM  VARIABLES 

THETA  » -A/ Al 
GC  TO  W 

COMPUTE  RATE  AVD  SCALING  PARAMETERS  FOR  HOMOGENEOUS 
PCISSON  PROCESS 
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3 PAR  = TINT*A 
I FLAG  = 1 
AINVRS  = l./A 

COMPUTE  NUMBER  OF  EXPONENTIAL  VARIATES  REQUIRED 

4 NMAX  * PAR+6.*SQRT (PAR ) 

IS  THIS  A HOMOGENEOUS  POISSON  PROCESS? 

IF  UFLAG.EQ.il  GO  TO  17 

GENERATE  REALIZATION  ON  POISSON  (THETA)  VARIATE 

5 CONTINUE 

CALL  PVAR  ( IS. THETAfM) 

IF  (M.EQ.O)  GO  TO  7 

CALCULATE  TIMES  OF  EVENTS 

CALL  SEXPON  (IS,T,NMAX) 

B * -Ai 

V = 0. 

JMAX  * NMAX+1 
OC  6 1*1, JMAX 

HAVE  NUMBER  OF  EVENTS  EXCEEDED  THE  MAXIMUM  NUMBER 
THAT  THE  ARRAY  CAN  HOLO? 

IF  (I.GT.NMAX)  GO  TO  8 

V = V+T ( I ) /((M-I+I )*B) 

IF  (V.GT.TINT)  GO  TO  9 
T ( I ) = V 

IF  U.EQ.M)  GO  TO  10 

6 CONTINUE 

NO  EVENTS  OCCURREO 

7 N = 0 
RETURN 

TOO  MANY  EVENTS  FOR  ARRAY.  INCREMENT  ERROR 
CODE  AND  TRY  AGAIN 

8 IER  = IER+l 
GO  TO  5 

THE  NUMBER  OF  EVENTS  OBSERVED  TO  OCCUR  IN  THIS 
NGN-HOMOGENEOUS  POISSON  PROCESS  IS  'N* 

S N = I-I 

IF  (N.EQ.O)  RETURN 
GO  TO  11 

10  N = M 

11  CONTINUE 

IS  THE  RATE  FUNCTION  INCREASING  OR  DECREASING? 

IF  (IFLAG.EQ.3I  GO  TO  13 

TIME  REVERSAL  TECHNIQUE  IS  NECESSARY 
DETERMINE  WHETHER  N IS  EVEN  OR  ODD 

ISIG  * MOD ( N, 2 ) 

NLOOP  = N/2 
N1  * N+i 

00  12  I *1 , NLOOP 
S = T(I) 

T(I)  * ER-T(Nl-I) 

T(N1-I)  = ER-S 
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C 


8 

C 


12  CONTINUE 

IF  (ISIG.EQ.l)  T( NL00P+1)=ER-T  (NLOOP+l) 

ARE  TIMES  OF  EVENTS  REQUESTED? 

IF  (II.EQ.O)  RETURN 
GO  TO  15 

IF  (II.NE.O)  GO  TO  15 
IF  (EL.EQ.O.)  RETURN 


13 

14 

15 

16 

17 

18 


DO  14  I «1 » N 
Till  « EL+Td) 

CONTINUE 
RETURN 

CALCULATE  TIMES  BETWEEN  EVENTS 
S * Till 

00  16  I *2t  N 
SI  * Till 
T ( I ) « T(  I ) -S 
S * SI 
CONTINUE 
RETURN 

THE  POISSON  PROCESS  IS  HOMOGENEOUS 

1 = 1 
U * 0. 

CALL  SEXPQN  ( I S t T tNMAX ) 

U = U+T( 1 1 

IF  (U.GT.PAR)  GO  TO  20 
I * 1+1 

IF  (I.GT.NMAX)  GO  TO  19 
GO  TO  18 

INCREMENT  ERROR  COOE 
IER  = IER+1 

TRY  AGAIN  WITH  NEW  STRING  OF  VARIATES 

GO  TO  17 
N = 1-1 

IF  (N.EQ.O)  RETURN 
IF  (II.EQ.l  ) GO  TO  22 

DO  21  I = 1 » N 

EL  = EL+AINVRS*TU) 

T(  I ) * EL 

21  CONTINUE 

RETURN 

22  00  23  I *1  » N 

T ( I ) « T( I ) *A INVRS 

23  CONTINUE 
RETURN 
ENO 


19 


20 


SUBROUTINE  PVAR  GENERATES  A POISSON  (THETA) 
VAR  I AT  Et  Mt  USING  THE  GAMMA  METHOD 

SUBROUTINE  PVAR  (IStTHETAtM) 

DIMENSION  T ( 5000 ) 

CCMMON  /HOLD/  T 


■i 

' 
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K * 0 
C * 16. 

D * .875 

IF  <THETA.LT.CJ  GO  TO  2 
GO  TO  5 
U = 1. 

CTN  * EXP(-THETAJ 

MMAX  * THET A+6 • *S QRT ( THET A) 

I * 1 

CALL  SRANO  (IS, T, MMAX) 

U « U*T ( I ) 

IF  (U.LT.CTN)  GO  TO  8 
I = 1 + 1 
K * K+l 

IF  (I.GT.MMAX)  GO  TO  3 
GO  TO  A 

NP  * INT(D*THETA) 

AN  * FLOAT ( NP ) 

CALL  GAMA  <AN,IS,G,1) 

IF  (G.GT. THETA)  GO  TO  6 

K * K+NP 

THETA  = THETA-G 

GO  TO  1 

U * THETA/G 

NP  = NP-1 

CALL  SRANO  <IS,T,NP) 


DO  7 1*1, NP 
IF  <T< I ) .LT.U) 
7 CONTINUE 


K*  1 


THE  VALUE  M IS  ASSUMED  BY  THE  POISSON  (THETA)  VARIATE 

8 M * K 
RETURN 
END 


1 


2 

3 


SUBROUTINE  REJECT  GENERATES  AN  ORDERED  SAMPLE 
OF  GIVEN  SIZE  FROM  THE  SECOND  COMPONENT 
OF  THE  ORIGINAL  INTENSITY  FUNCTION 
USING  A REJECTION-ACCEPTANCE  ALGORITHM 

SUBROUTINE  REJECT  ( I S , EL ,CSTAR , PVEC ,Q ,L ) 
DIMENSION  V (500) , PVEC(5) 

01 MENSI ON  T ( 5000) 

CCMMON  /HOLD/  T 
L20  = L*1 0 

IF  (L20.GT.500)  L20»500 
LI  » L+l 
K * 1 
J = 0 

CALL  SR  AND  (IS,V,L20) 

DO  2 I * 1,  L20 
J * J+l 

T(K)  « Q*V( J) +EL 
J * J+l 

IF  (V(J).LT.CALC(PVEC,T(K))/CSTAR)  K*K* 1 
IF  (K.EQ.Ll)  GO  TO  3 
IF  (J.GE.L20-1)  GO  TO  I 
CONTINUE 

IF  (K.LT.L)  GO  TO  1 
CALL  PXSORT  <T,1,L) 

RETURN 

ENO 
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SUBROUTINE  HELP  EVALUATES  THE  INTEGRATED 
FUNCTION  OVER  THE  INTERVAL  ( EL « ER ) 

SUBROUTINE  HELP  ( A , Al , A2, EL ,ER , XX) 

DOUBLE  PRECISION  MMDAW » BB » A A 
Z = SORT ( ABSCA2) ) 

Y » (Al*Z)/(2.*A2 ) 

AA  » Z*  EL+Y 
BB  = Z*ER+Y 
CC  * A-A1*A1/(4.*A2) 

CC  » EXP< CC )/Z 

IF  (A2.LT.O.)  30  TO  l 

Q I » DEXP(AA**2)*MMDAW(AA) 

Q2  = DEXP( BB**2 )* MMDAW ( BB) 

XX  » CC*(Q2H31) 

RETURN 

CC  = CC*. 8862269 

XX  = CC*(OERF(BB ) -OER  F{ AA ) ) 

RETURN 

END 


INTENSITY 


SUBROUTINE  COLATE  SUPERPOSES  TWO  ORDERED 
EVENT  STREAMS  OVER  THE  SAME  INTERVAL 

SUBROUTINE  COLATE  ( LBGN . LF I N, JBGN) 
DIMENSION  T IMES  ( 5000 ) t T(5000) 

COMMON  /MIKE/  TIMES/HOLD/T 
I a JBGN 
J a I 
K = LBGN 

1 IF  (TIMES(I).LT.TIMES(K) ) GO  TO  2 
T ( J ) » TIMES(K) 

J = J+l 
K a K + l 

IF  (K.GT.LFIN)  GO  TO  3 
GO  TO  1 

2 T ( J)  a TIMES(I) 

J = J + l 

I » I + l 

IF  ( I. EQ.LBGN)  GO  TO  5 
GO  TO  l 

3 II  = LBGN-1 

DO  4 Nal.II 
T ( J ) » TIMES(N) 

J « J+l 

4 CONTINUE 
RETURN 

5 CONTINUE 

00  6 N*K»  LF  IN 
T ( N ) a TIMES(N) 

6 CONTINUE 

RETURN 

ENO 


FUNCTION  CALC  EVALUATES  THE  SECOND  COMPONENT  OF 
THE  DECOMPOSED  INTENSITY  FUNCTION 
FOR  ANY  INPUT  VALUE. 

FUNCTION  CALC  (P,ABSA) 

DIMENSION  P ( 5) 

X » P(l)+P(2)*ABS  A+P ( 3 )*A3S A**2 
XX  * P<  4) +P ( 5 ) * AB SA 
CALC  » EXP ( X) “EXP (XX ) 

RETURN 

ENO 
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